The mechanistic model of neurovascular coupling was developed and studied by Sten et al. (2020). This model describes and predicts the arteriolar dilation data of mice under various stimulations while anaesthetised and awake. We reconstructed the model in CellML, using a modular approach for each neuronal pathway, and successfully reproduced the original experiments (see the figures in this article and Sten et al. (2020)). With the success of the result reproduction, the CellML model can now be injected into other OpenCOR workflows to obtain a mechanistic hemodynamic response function of neurovascular coupling arteriolar dilation.
Neurovascular coupling (NVC) is the cerebral pathway to modulate blood flow by manipulating local vascular tone Iadecola, 2017Schaeffer & Iadecola, 2021. This modulation is required based on the metabolic demand of the tissue, and its impairment can lead to neurodegeneration and cognitive disease Sweeney et al., 2019Iadecola & Gottesman, 2019. To prevent this impairment, we need a detailed understanding of the mechanisms of the NVC unit. In line with this goal, researchers have developed mechanistic mathematical models of NVC Hart et al., 2019Sten et al., 2020 based on experiment. The mechanistic hemodynamic response functions (HRF) of the models under various stimulations are then tested against in vivo experiments. With success, the model is then used to explore other core predictions produced increasing our knowledge of NVC and motivating further study. The model presented in Sten et al. (2020) considers GABAergic interneurons for the first time and is capable of accurately modelling arteriolar dilation data from mice under awake and anaesthetised conditions, as reported by Uhlirova et al. (2016).
Here, we present a CellML version of the model using the system of equations described in Sten et al. (2020) without modifications. Our CellML implementation is modular, allowing future studies to import separate pyramidal (Pyr), GABAergic Nitric Oxide (NO) and GABAergic Neuropeptide Y (NPY) control. With this model, our objective is to reproduce all steady-state variables that optimise arteriolar dilation results in Sten et al. (2020). We omit the reproduction of confidence intervals generated by Markov Chain Monte Carlo (MCMC) simulations, since the goal is model reproducibility and reusability, not evaluation of model fidelity or sensitivity, as provided in the primary publication Sten et al., 2020.
The proposed CellML implementation can be run in OpenCOR and Python software, extending the implementation of the primary publication in MATLAB R2017b. This implementation enables more researchers to access the model.
The model of Sten et al. (2020) is a mechanistic model of neurovascular coupling that begins with neuronal stimulation and ends with modulation of arterial tone. The model was designed based on the iterative philosophy of systems biology, which allows quantitative experimental data to drive model development and improvement. If a proposed model can capture experimental data and is physiologically sound, then the model is accepted.
At times, an accepted model may be cumbersome due to the number of parameters and the complexity of the system. In these instances, complexity can be reduced if the model is minimised. This is the process of removing parts of the model and re-optimising to the data. This continues until the minimised model can no longer capture the data with confidence Lundengård et al., 2016. The model presented in the primary publication is a minimised systems biology model. Interpretations of state variables can be made, but it is important to recognise that potentially multiple processes with similar responses may be lumped into one set of equations.
The model is divided into three sections, neuronal, intracellular, and vascular. In the neuronal section, each type of neuron (pyramidal, GABAergic NO, and GABAergic NPY) can receive stimulation that increases their activity. This activity can activate the other two neurons in the case of pyramidal activation or inhibit activity in the case of GABAergic interneurons. In the respective signalling pathways of neurons, activity stimulates depolarisation and a subsequent increase in calcium that activates the signalling pathway to produce vasodilatory (pyramidal and GABAergic NO pathways) or vasoconstrictive (GABAergic NPY path) factors. These final three outputs are collected and summed with respective weighting coefficients to simulate the response of the vessel. The readers are referred to the primary publication for complete details of the model Sten et al., 2020.
The model contains a total of coupled ODEs based on mass-action, or Michaelis-Menten kinetics. The original model was implemented in MATLAB 2017b and used the third-party package AMICI to run the CVODES solver. The model was optimised using MCMC methods provided by the MEIGO and PESTO packages. The scripts to run the original model are provided in the supplementary material of Sten et al. (2020).
The structure of the model is built following recommended CellML best practises (see Figure 1). This includes modularisation of the pyramidal, GABAergic NO and GABAergic NPY neuronal signalling pathways for best combined prediction or separate use. Best practises also encourage the separation of model mathematics and parameters so that one model can handle various parameter sets. For this model, the parameters simulate different neuronal dynamics. Examples that will impact neuronal dynamics include the type of stimulation, the conditions of the experiment (anaesthesia), and the genetic modulation that manifests itself as coefficient modification. We have achieved these criteria, as well as reduced the number of input parameters by including a toggle for anaesthesia versus awake conditions in our parameter files. This is a benefit compared to the five separate models required in the primary publication. This simplification also reduces the toll of scripts, allowing users to run different scenarios by only altering the set of parameters. To reproduce all results, only six parameters in the parameters file need to be modified. These six parameters define the start and stop times for stimulation, whether the unit is “awake” or “anaesthetised”, and finally the stimulation scaling constants are discussed in detail in the primary publication.
The CellML model was executed with the optimised rate constants and other parameters tabulated in the supplementary material of the primary publication. These values are rounded from what was optimised in the original code; differences if any, will be discussed in the results. The CellML model is initialised using steady-state values from its own simulation run for seconds and not from the tabulated in the supplementary material. This accounts for any steady-state result errors that may have been caused by rounding. Again, the differences will be discussed in the results.
This model will be used to reproduce Figures , and of the primary publication.
2.3Clarifications and Modifications¶
In Figure of the primary publication, the stimulation windows are shown in each graph, but the exact timing of neuronal stimulation is not presented in the manuscript or the supplementary materials. It is presented in the code, which we tabulate here for clarity.
Table 1:Stimulation times to reproduce figures.
|Figure #||Model Conditions||Stimulation||Stimulation|
|(Primary)||Start Time (sec)||Length (sec)|
|Figure 3||Awake-Sensory Stim||0||1|
|Figure 5/6||Awake-Optogenetic (GABAergic)||0||0.4|
In the reproduction of Figure c, the label in is ; however, the normalisation to reproduce this plot followed Equation 1 below, which means .
Finally, to reproduce the plot data in Figure of the primary publication, the first row of the plots corresponding to the neuronal activity of the NPY and NO pathways is not normalised as suggested in the primary publication. We have changed the title to “Cell Activity”, which is unitless.
The results were generated using the 2022-05-31 snapshot of OpenCOR Garny & Hunter, 2015 and MATLAB R2021a to produce the plots . The MATLAB code requires software newer than R2019b. The entire code is available at https://
Table 2:Parameter set names
|Start Text||Simulation Tag (%%)||End Text|
3.1Step1: Load parameters¶
Open sten2020_NVC_main.cellml and on line alter the parameter name to any of the tags (%%) in Table 2. All parameter filenames are in the same location as the main function and follow the naming of the figure to be reproduced from the primary publication.
3.2Step2: Simulation Parameters¶
The equations were solved using the CVODE solver with a maximum timestep-size equal to . The model is already initialised in steady state. All run times are seconds with a timestep of seconds.
When the simulation is complete, export all the data to .csv. Use the default name sten2020_NVC_main_data.csv but amend the parameter case (%%) to the end as shown sten2020_NVC_main_data_
%%.csv. Place each csv file into an <Output CSVs> folder in the workspace. This is filled with the results by default.
3.4Step4: Generate Plots¶
To reproduce the figures, run the chosen Matlab function in the folder <MatlabFiles> from within its directory. The name of the script corresponds to the figure it produces; for example, running
will produce Figure in the primary publication.
Starting with steady-state values, only values were not identical to those published in the supplementary material of the primary publication. The steady-state values for arachidonic acid, NO, and NO were , , and in the CellML model and , , and in the primary publication; all differences are less than %. This led us to conclude that using the rounded numbers was safe to continue with figure reproduction.
The first figure reproduced is Figure of the primary publication shown in Figure 2. In this figure, the results of arterial dilation are shown for different stimulus cases. The parameters used to produce this graph come from the best-fit parameters in the primary publication. All subplots agree with the primary publication.
The second figure reproduced is Figure of the primary publication shown in Figure 3. This figure emphasises the impact of using Michaelis-Menten kinetics for the production of NPY in vascular smooth muscle (vsm). All subplots agree with the primary publication.
The last reproduced figure is Figure of the primary publication shown in Figure 4. This figure shows the difference in all state variables of the GABAergic NO and NPY pathways for awake and anaesthetised conditions, highlighting that the prolonged undershoot of the anaesthetic response is due to the Michalis-Menten kinetics of NPY. It also shows the impact of NO and its initial vasodilatory effect. For a complete interpretation and all acronyms, see the primary publication. All sub-plots agree with the primary publication.
In this paper, we present the CellML version of an NVC model originally developed by Sten et al. (2020) for NVC. Our implementation is modular and can be tuned to many coupling scenarios without the main code modification. Figures illustrating the model predictions from the primary publication were produced without error, and only minor clarifications are made to the titles in the primary publication’s figures.
- Sten, S., Elinder, F., Cedersund, G., & Engström, M. (2020). A quantitative analysis of cell-specific contributions and the role of anesthetics to the neurovascular coupling. NeuroImage, 215. 10.1016/j.neuroimage.2020.116827
- Iadecola, C. (2017). The Neurovascular Unit Coming of Age: A Journey through Neurovascular Coupling in Health and Disease. Neuron, 96(1), 17–42. 10.1016/j.neuron.2017.07.030
- Schaeffer, S., & Iadecola, C. (2021). Revisiting the neurovascular unit. Nature Neuroscience, 24(9), 1198–1209. 10.1038/s41593-021-00904-7
- Sweeney, M. D., Montagne, A., Sagare, A. P., Nation, D. A., Schneider, L. S., Chui, H. C., Harrington, M. G., Pa, J., Law, M., Wang, D. J. J., Jacobs, R. E., Doubal, F. N., Ramirez, J., Black, S. E., Nedergaard, M., Benveniste, H., Dichgans, M., Iadecola, C., Love, S., … Zlokovic, B. V. (2019). Vascular dysfunction—The disregarded partner of Alzheimer’s disease. Alzheimer’s& Dementia, 15(1), 158–167. 10.1016/j.jalz.2018.07.222
- Iadecola, C., & Gottesman, R. F. (2019). Neurovascular and Cognitive Dysfunction in Hypertension. Circulation Research, 124(7), 1025–1044. 10.1161/CIRCRESAHA.118.313260