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

1Introduction¶

Lees-Green et al. (2014) 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 & Hunter, 2015 Python interpreter to perform the simulations and analyses. In this manuscript, we focus on reproducibility and reusability.

MATLAB Shampine & 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 Section 3.2.

2Model description¶

2.1Primary 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 (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. (6) is used to calculate the cytosolic Ca$^{2+}$ concentrations for all the above mentioned model variations.

The High-Cl(CaV: Voltage-Activated Ca$^{2+}$ channel) model variation incorporates a long-lasting voltage-gated Ca$^{2+}$ channel (I$_{CaV}$), so the total ionic current across the membrane is following (7). We used (8) to determine the change in cytosolic Ca$^{2+}$.

In the case of the Low-Cl(NSCC: Ca$^{2+}$ Activated Non-Selective Channel), (9) and (10) are used to calculate the total ionic current and the change in cytosolic Ca$^{2+}$, respectively.

The model constituents vary according to the different model configurations, in regard to chloride concentrations: high chloride environment ($E_{Cl} = - 20.2$ mV, $C_{Cl} = 78$ mM), and low chloride environment ($E_{Cl} = - 49.7$ mV, $C_{Cl} = 25.85$ mM). See Table 3 for details on the model configurations.

2.2Model 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).

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 ( Fig4-1a & Fig4-1b) to hold the data. For more detailed information, see Section 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.

3Model Modifications¶

3.1Mathematical Equations¶

The equations are the same as reported in the supplementary material of Lees-Green et al. (2014) with the exception of some of the flux definitions. There are no definitions of the fluxes for some conductances: voltage-gated calcium channel ($J_{CaV}$), T-type calcium channel ($J_{CaT}$), Soc channel ($J_{Soc})$, and voltage-gated nonselective channels ($J_{NSCC}$). Definitions are defined (see (11)) from the references and compared to the similar descriptions of the other fluxes with similar behaviour in the original work and confirmed through the MATLAB updated model of ICC. The following equation defines the total fluxes for the conductances mentioned above:

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

3.2Parameter 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 ($\tau_{Ano1}$) was not defined in the primary paper; we highlight that the rescaling factor was needed to produce Figure 3 from Lees-Green et al. (2014).

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 Section 3.3.1.

3. To reproduce Figure 5A-5C, the $K$ parameter for voltage-dependent gating equation ($F_{NaV}$) 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.

4. To reproduce 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 Section 3.3.2.

3.3Computational Simulation¶

The differential$/$algebraic equations system is solved using solver CVODE in CellML. with the relative tolerance and time step of $Retol=e-7$ and $step =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).

3.3.1Ano1 model¶

The Ano1 model is defined to study the kinetics of the Ano1 channels (g = 0.5 nS, E$_{Cl}$ = 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 $\mu M$ (Figure 3A), 1 $\mu M$ (Figure 3B), and 10 $\mu M$ (Figure 3C), see (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 $\mu$ M to 30 $\mu$ M for 800 ms. For more detailed information see (13) and Table 2.

3.3.2ICC model¶

The ICC model not only studies the effect of altering E$_{Cl}$ 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.

The high chloride environment contains four different variations, as the following:

1. HIGH-CL(NSV): Non-Selective Voltage Activated Channel, to reproduce Figure 4B and F; ‘HCl-NSV.json’ is loaded to update the model’s parameters and initial conditions accordingly. Data is saved in the CVS files: Fig4-2a and Fig4-2b. All the data held in the appropriate string of formats ‘ia’ and ‘ib’ where ‘i’ indicates the figure subplot, ‘a’ and ‘b’ indicate the data related to the Wild-type, and Ano1 knockout scenarios, respectively.

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.

3. High-Cl(NSCa): Ion currents specific to Ca2-activated non-selective channel, to reproduce Figure 4C and G; ‘HCl-NSCC.json’ is loaded.

4. High-Cl(CaV): Ion current specific to Voltage-gated Ca channel, to reproduce Figure 4D and H; ‘HCl-CaV.json’ is loaded.

The low chloride environment consists of two different variations:

1. Low-Cl(NaV): Ion currents specific to Voltage-gated Na channel, to reproduce 5A-5C; ‘LCl-NaV.json’ is loaded.

2. Low-Cl(NSCa): Ion currents specific to Ca2-activated non-selective channel, to reproduce Figure 5B and 5D; ‘LCl-NSCC.json’ is loaded.

3.4Model Results¶

We employed the Engauge Digitizer software Mitchell et al., 2020 to extract the data from all figures from Lees-Green et al. (2014) to present alongside the results of the present work. The reproduction of all figures of Lees-Green et al. (2014) is given in Figures 1-3, which show good agreement against the data of the primary paper. Solid lines show the output from this curation effort, and stars show discrete points found by the Engauge Digitizer of the figures in the primary paper.

4Discussion¶

We demonstrated that most results of Lees-Green et al. (2014) can be reproduced given some added equations (Section 3.1) and adjustments (Table 1) made to model parameter values. We implemented the model Lees-Green et al., 2014 using CellML, replicated all the model behaviors reported in the primary paper, and summarized detailed experiment settings for simulation in Section 3.3. While we are generally able to reproduce the simulation results there are some cases where the results presented here are not identical. The reason for this is the lack of information about the initial conditions from the primary paper leading to slight phase shifts in the oscillatory behaviour as the simulations are not starting from the same initial state. These differences do not impact their interpretation and the conclusions drawn in the primary 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.