Reproducibility Study for a Computational Model of the Neurovascular Coupling Unit

Abstract¶

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.

Keywords:CellMLneurovascular couplinganaesthesiareproducibility¶

1Introduction¶

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^[1] and Python software, extending the implementation of the primary publication in MATLAB R2017b. This implementation enables more researchers to access the model.

2Model Description¶

2.1Primary Publication¶

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 $13$ 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^[2] to run the CVODES solver. The model was optimised using MCMC methods provided by the MEIGO^[3] and PESTO^[4] packages. The scripts to run the original model are provided in the supplementary material of Sten et al. (2020).

2.2CellML Model¶

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 $500$ 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 $3$ , $5$ and $6$ of the primary publication.

Figure 1:Arrangement of the model showing the CellML scripts involved. Sten2020_NVC_main.cellml imports components from the other scripts and constructs the model through mapping of variables between components.

2.3Clarifications and Modifications¶

In Figure $3$ 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
	Awake-Optogenetic (GABAergic)	0	0.4
	Anesthesia-Sensory	0	2
	Anesthesia-Optogenetic (GABAergic)	0.55	0.45
	Anesthesia-Optogenetic (Pyramidal)	0.9	0.1
Figure 5/6	Awake-Optogenetic (GABAergic)	0	0.4
	Anesthesia-Optogenetic (GABAergic)	0	0.45

In the reproduction of Figure $5$ c, the label in $\Delta v_1=1$ is $V_\mathrm{max}$ ; however, the normalisation to reproduce this plot followed Equation 1 below, which means $1\equiv V_\mathrm{max}-v_1^\mathrm{SteadyState}$ .

\Delta v_1 = \frac{\frac{V_\mathrm{max}NPY(t)}{K_{M}+NPY(t)}-v_1^{\mathrm{SteadyState}}}{V_\mathrm{max}-v_1^{\mathrm{SteadyState}}}

(1)

v_1^\mathrm{SteadyState} = \frac{V_\mathrm{max}NPY_\mathrm{SteadyState}}{K_{M}+NPY_\mathrm{SteadyState}}

(2)

Finally, to reproduce the plot data in Figure $6$ 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.

3Model Execution¶

The results were generated using the 2022-05-31 snapshot of OpenCOR Garny & Hunter, 2015 and MATLAB R2021a to produce the plots The Mathworks, Inc., 2021. The MATLAB code requires software newer than R2019b. The entire code is available at https://models.physiomeproject.org/workspace/8a2. To generate results, follow this step-wise procedure for each parameter set named in Table 2.

Table 2:Parameter set names

Start Text	Simulation Tag (%%)	End Text
	Fig3A	.cellml
	Fig3B	.cellml
	Fig3C	.cellml
sten2020_parameters_	Fig3D	.cellml
	Fig3E	.cellml
	Fig5and6_Ana	.cellml
	Fig5and6_Awake	.cellml

3.1Step1: Load parameters¶

Open sten2020_NVC_main.cellml and on line $24$ 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 $0.01$ . The model is already initialised in steady state. All run times are $50$ seconds with a timestep of $0.01$ seconds.

3.3Step3: Export¶

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 Figure3.m will produce Figure $3$ in the primary publication.

4Model Results¶

Starting with steady-state values, only $3$ 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 $_\mathrm{vsm}$ were $6.224$ , $2944.421$ , and $352.127$ in the CellML model and $6.226$ , $2944.636$ , and $352.139$ in the primary publication; all differences are less than $0.01$ %. This led us to conclude that using the rounded numbers was safe to continue with figure reproduction.

The first figure reproduced is Figure $3$ of the primary publication shown in Figure 2. In this figure, the results of arterial dilation are shown for $5$ 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.

Best estimated model simulation of arteriolar dilation to sensory or optogenetic stimulation in awake and anesthetised mice. This figure matches the subfigures B-D-F-H-J of Figure 3 in the primary publication and can be reproduced using Figure3.m file. — Figure 2:Best estimated model simulation of arteriolar dilation to sensory or optogenetic stimulation in awake and anesthetised mice. This figure matches the subfigures B-D-F-H-J of Figure $3$ in the primary publication and can be reproduced using `Figure3.m` file.

The second figure reproduced is Figure $5$ 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.

NPY response from optogenetic stimulation of GABAergic interneurons in awake and anesthetised mice. (top left) Change in NPY from baseline steady state for awake and anesthetised optogenetic GABAergic stimulation. (top right) Change in rate of NPY_\mathrm{vsm} production as a function of NPY produced from baseline showing minor changes in rate for large increases in production in anesthetised vs awake conditions. (bottom) Change in NPY_\mathrm{vsm} over time for anesthetised and awake stimulation. This figure matches the subfigures B-C-D of Figure 5 in the primary publication and can be reproduced using Figure5.m file. — Figure 3:NPY response from optogenetic stimulation of GABAergic interneurons in awake and anesthetised mice. (top left) Change in NPY from baseline steady state for awake and anesthetised optogenetic GABAergic stimulation. (top right) Change in rate of NPY $_\mathrm{vsm}$ production as a function of NPY produced from baseline showing minor changes in rate for large increases in production in anesthetised vs awake conditions. (bottom) Change in NPY $_\mathrm{vsm}$ over time for anesthetised and awake stimulation. This figure matches the subfigures B-C-D of Figure $5$ in the primary publication and can be reproduced using `Figure5.m` file.

The last reproduced figure is Figure $6$ 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 $_\mathrm{vsm}$ . 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.

Complete presentation of the NO and NPY state variables for optogenetic stimulation of GABAergic interneurons in awake (red) and anesthetised (black) mice. This figure matches Figure 6 in the primary publication and can be reproduced using Figure6.m file. — Figure 4:Complete presentation of the NO and NPY state variables for optogenetic stimulation of GABAergic interneurons in awake (red) and anesthetised (black) mice. This figure matches Figure $6$ in the primary publication and can be reproduced using `Figure6.m` file.

5Discussion¶

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.

References¶

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