OPEN ACCESS Reproducible Model

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 ﬁgures 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 workﬂows to obtain a mechanistic hemodynamic response function of neurovascular coupling arteriolar dilation.


Introduction
Neurovascular coupling (NVC) is the cerebral pathway to modulate blood flow by manipulating local vascular tone (Iadecola, 2017;Schaeffer and 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., 2019;Iadecola and 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., 2019;Sten 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 OpenCOR1 and Python software, extending the implementation of the primary publication in MATLAB R2017b.This implementation enables more researchers to access the model.

Primary 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 AMICI2 to run the CVODES solver.The model was optimised using MCMC methods provided by the MEIGO3 and PESTO4 packages.The scripts to run the original model are provided in the supplementary material of Sten et al. (2020).

CellML 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.

Clarifications 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.
In the reproduction of Figure 5c, the label in ∆v 1 = 1 is V max ; however, the normalisation to reproduce this plot followed Equation 1 below, which means 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 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.

Model Execution
The results were generated using the 2022-05-31 snapshot of OpenCOR (Garny and Hunter, 2015) and MATLAB R2021a to produce the plots (Mat, 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. .cellmlFig5and6_Awake .cellml

Step1: 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.

Step2: 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.

Step3: 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. 4/8

Step4: 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.

Model 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 vsm were 6. 224, 2944.421, and352.127 in the CellML model and6.226, 2944.636, and352.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.is due to the Michalis-Menten kinetics of NPY 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.

Discussion
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.

Figure 1 .
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.

Figure 2 . 8 Figure 3 .
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 ofFigure3in the primary publication and can be reproduced using Figure3.m file.

Figure 4 .
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 matchesFigure6in the primary publication and can be reproduced using Figure6.m file.

Table 2 .
Parameter set names