Computational modeling of anoctamin 1 calcium-activated chloride channels as pacemaker channels in interstitial cells of Cajal

The system of equations and simulation results presented in Lees-Green et al. (2014) are veriﬁed and reproduced in the current paper. To demonstrate reproducibility, we here describe the model encoded in CellML and document the diﬀerences between our curated model and that published by Lees-Green et al..


Introduction
describes a biophysical computational model of anoctamin 1 (Ano 1) calcium-activated chloride channels (pacemaker channels in interstitial cells of Cajal). The current work involves the mathematical curation of the model in the CellML format (Cuellar et al., 2003), which can reproduce the behaviour of each of the model's variations presented in Lees-Green et al. (2014). Figure 2 in Lees-Green et al. (2014) is not included as it does not represent the introduced model's results. Additionally, some inconsistencies have been identified and are discussed in this work. A persistent workspace for this work is available in the Physiome Model Repository at https://models.physiomeproject.org/workspace/83c. The specific version of this implementation used to produce the presented results is included in the OMEX Archive (Bergmann et al., 2014).
This implementation includes the Python source code to generate the model simulations in a 'Simulations' directory. The corresponding CellML implementation is available in a 'Experiments' directory containing the encoded mathematical model. This CellML file on its own does not reproduce the figures due to limited solver capabilities at time of testing. As such, we have elected to describe the required simulation experiments in Python and use the OpenCOR (Garny and Hunter, 2015) Python interpreter to perform the simulations and analyses. In this manuscript, we focus on reproducibility and reusability.
MATLAB (Shampine and Reichelt, 1997) code was obtained from the author of Lees-Green et al. (2014) to reproduce the figures of that paper. While there was no exact correspondence between equations presented in the original article and the provided MATLAB code, all the figures (corresponding to the model results) were produced, and numeric values were matched with some modifications to a few parameter values; for more details, please see subsection 3.2.

Primary Publication
The model uses a Hodgkin-Huxley type formulation. The cell membrane lipid bilayer is represented as a capacitance (C m ), and the ion channels in the membrane are represented as conductance. The change in the transmembrane potential (V m ) over time depends on the sum of the individual ion currents through each class of ion channel in the cell, see Equation 1, There are 12 different ion conductances in this model (to see the components of the model, see Table 3 and the equations below). We noted that the total ionic current and total flux vary according to the different model variations. For the following model variations: High-Cl(NaV: voltage activated Na Channel), High-Cl(NSV: Voltage Activated Non-Selective Channel), High-Cl(NSCC: Ca 2+ Activated Non-Selective Channel), and Low-Cl(NaV: Voltage Activated Na + Channel); Equations 2, 3, 4, and 5 represent the total ionic current, respectively. Equation 6 is used to calculate the cytosolic Ca 2+ concentrations for all the above mentioned model variations.
I t ot (HCl-NaV) = I SO C + I Ano1 + I C aT + I BK + I B N a + I N aV , I t ot (HCl-NSCC) = I SO C + I Ano1 + I C aT + I BK + I B N a + I N S C C , I t ot (LCl-NaV) = I SO C + I Ano1 + I C aT + I BK + I B N a + I N aV + I KV + I K E RG , The High-Cl(CaV: Voltage-Activated Ca 2+ channel) model variation incorporates a long-lasting voltage-gated Ca 2+ channel (I C aV ), so the total ionic current across the membrane is following Equation 7. We used Equation 8 to determine the change in cytosolic Ca 2+ .
I t ot (HCl-CaV) = I SO C + I Ano1 + I C aT + I BK + I B N a + I C aV , In the case of the Low-Cl(NSCC: Ca 2+ Activated Non-Selective Channel), Equation 9 and Equation 10 are used to calculate the total ionic current and the change in cytosolic Ca 2+ , respectively.
I t ot (LCl-NSCC) = I SO C + I Ano1 + I C aT + I BK + I B N a + I N S C C + I N C X , The model constituents vary according to the different model configurations, in regard to chloride concentrations: high chloride environment (E C l = −20.2 mV, C C l = 78 mM), and low chloride environment (E C l = −49.7 mV, C C l = 25.85 mM). See Table 3 for details on the model configurations.

Model Implementation
The implementation of the model was performed using CellML (Cuellar et al., 2003). Some simulation values stated in Lees-Green et al. (2014) resulted in different model behaviour. We adjusted them as outlined in Table 1. The simulation experiments presented here were performed using OpenCOR (Snapshot 2021-10-5).

2/9
The simulation experiments in Lees-Green et al. (2014, Figures 4 & 5 ) were produced in various stages corresponding to the model variations, as summarised in Table 3. We implemented the model configuration for each variation using the related parameters input file for each. For example, to produce the data related to the model variation 'H-Cl (NaV)', we load the model configuration file 'HCl-NSV.json' and we choose the CSV files  to hold the data. For more detailed information, see subsubsection 3.3.2 and the simulation folder in the mentioned workspace ('ICC-Lees-Green.py'). Then, we executed the Python script to plot the related data corresponding to each figure. Experimental conditions specific to each variations are given in Table 3.

Mathematical Equations
The where I and V indicate the current and the cell volume, respectively.

Parameter Values
The baseline parameter values are provided in the supplementary material of Lees-Green et al. (2014). In an attempt to reproduce the same results, we applied a few modifications as listed below (also see Table 1): 1. The rescaling factor of time constant for the Ano1 channel (τ Ano1 ) was not defined in the primary paper; we highlight that the rescaling factor was needed to produce 2. In Figure 3D, Lees-Green et al. defined the step changes in calcium concentrations. Applying the mentioned range in the primary paper could not reproduce the same result; for more information, see subsubsection 3.3.1.
3. To reproduce Figure 5A-5C, the K parameter for voltage-dependent gating equation (F N aV ) is changed (from -4.5 mV to 4.5 mV). The K parameters corresponding to the inactivation of the voltage-gated channels are always positive. Figure 5B-5D, the conductance value for the 'CaV' channel is needed to be reduced from the original value g = 4 nS to g = 3.72 nS.

5.
No initial conditions were defined in the primary paper. All initial conditions were extracted from the original article using the Engauge Digitizer software (Mitchell et al., 2020). These values are embedded in the 'json' files regarding each variation; see subsection subsubsection 3.3.2.

Computational Simulation
The differential/algebraic equations system is solved using solver CVODE in CellML. with the relative tolerance and time step of R et ol = e − 7 and st ep = 1e − 3, respectively. In this work, the model results are categorized into two different parts: The Ano1 model (Figure 1) and the ICC model ( Figure 2 and Figure 3).

Ano1 model
The Ano1 model is defined to study the kinetics of the Ano1 channels (g = 0.5 nS, E C l = 0 mV). In the voltage-clamp simulations, the voltage protocol stepped from a holding potential of 0 mV to a membrane potential (V m ) ranging between -100 and +100 mV in 40 mV increments for 800 ms, followed by a hyperpolarizing step -100 mV. Voltage clamp simulations were carried out at three different values of fixed [Ca 2+ ]: 0.1 µM ( Figure 3A), 1 µM ( Figure 3B), and 10 µM ( Figure 3C), see Equation 12.
V t est 1s < t <= 9s −100 mV ot her w i se (12) Figure 3D was reproduced using the Ano1 model by setting the membrane potential (V m ) at a fixed holding potential (ranging from -100 to +100 mV in 40 mV increments) and applying a step-change in Ca Ano1 from 15 µ M to 30 µ M for 800 ms. For more detailed information see Equation 13 and Table 2.
C a Ano1 = C a t est 1s < t i me <= 9s 0 µM ot her w i se

ICC model
The ICC model not only studies the effect of altering E C l but also investigates and tests whether a variety of different ion currents are plausible candidates for the plateau current that maintains the plateau phase of the slow-wave. Simulations are divided into two distinct categories: high chloride and low chloride environments. 2. High-Cl(NaV): Ion current specific to Voltage-gated Na channel, to reproduce Figure 4A and E; the JSON file 'HCl-NaV.json' is loaded.

Discussion
We demonstrated that most results of Lees-Green et al. (2014) can be reproduced given some added equations (subsection 3.1) and adjustments (Table 1) Figure 4), use the Python script 'ICC _ Lees _ Green.py' to produce the simulation results and 'Plot _ Fig4 _ ICC.py' to plot the results. A-H: membrane potentials were simulated using the wild-type (WT) and Ano1 knockout (KO) scenarios (blue and red lines, respectively).
paper of Lees-Green et al. (2014). Python code was needed in some cases and can be found at https://models.physiomeproject.org/workspace/762. As this paper follows the tenets of FAIR (Wilkinson et al., 2016), further credibility is lent to the primary models.

Box 1: Criteria for repeatability and reproducibility
Model source code provided: Source code: a standard procedural language is used (e.g. MATLAB, Python, C) There are details/documentation on how the source code was compiled There are details on how to run the code in the provided documentation The initial conditions are provided for each of the simulations Details for creating reported graphical results from the simulation results Source code: a declarative language is used (e.g. SBML, CellML, NeuroML) The algorithms used are defined or cited in previous articles The algorithm parameters are defined Post-processing of the results are described in sufficient detail

Executable model provided:
The model is executable without source (e.g. desktop application, compiled code, online service) There are sufficient details to repeat the required simulation experiments

The model is described mathematically in the article(s):
Equations representing the biological system There are tables or lists of parameter values There are tables or lists of initial conditions

Machine-readable tables of initial conditions
The simulation experiments using the model are described mathematically in the article: Integration algorithms used are defined Stochastic algorithms used are defined Random number generator algorithms used are defined Parameter fitting algorithms are defined The paper indicates how the algorithms yield the desired output