Abstract

The system of equations and figures presented in Saucerman et al. (2003) are verified and reproduced in this paper’s curation effort. These checks are performed in MATLAB. With some parameter additions and modifications, we can reproduce all figures with small mismatches.

Keywords:beta-adrenergiccardiac myocyteMATLAB

1Introduction

β\beta-adrenergic signalling within the cardiac myocyte is a well-documented signalling pathway, but few mathematical models exist of the complete system. Saucerman et al. (2003) published such a model of the activation of the signalling cascade from a ligand binding the receptor on the cell surface, to triggering the activation of the beta-1 adrenergic receptor (β\beta1AR), cyclic AMP (cAMP), adenylyl cyclase (AC), protein kinase A (PKA), and phospholamban (PLB). Subsequent downstream effects on ion channel activity are the result, notably through modifying the activity of calcium channels in the regulation of an action potential.

This present work involves the mathematical curation the model. MATLAB MATLAB, 2017 code was obtained from the first author of Saucerman et al. (2003) in an effort to reproduce the figures of that paper. While there is not an exact correspondence between equations presented in the paper and code provided, the figure trends are still reproduced, and numeric values are matched with some modifications/additions made to parameter values, as detailed in under Model Modifications.

A persistent workspace of the code is available in the Physiome Model Repository, at https://models.physiomeproject.org/workspace/6fa.

This workspace includes the MATLAB source code for running the model. The corresponding CellML Cuellar et al., 2003 implementation is available at https://models.physiomeproject.org/e/6f6 which contains the encoding of the mathematics, but does not reproduce the figures due to limited solver capabilities at time of testing. As such, we have elected to do this present curation in MATLAB.

2Model description

The model is broken into several modules:

  • β\beta1AR activation

  • Gs protein activation

  • cAMP metabolism

  • PKA activation

  • PLB regulation

  • L-type calcium channel (LCC) regulation

Once the cell is activated by an agonist isoproterenol (Iso), the collective effects of these modules are studied downstream. This includes changes to ion channel activity towards the production and management of an action potential.

3Model Modifications

3.1Mathematical Equations

These are the same as reported in the supplementary material of Saucerman et al. (2003) with the exception of the L-type calcium channel. This is a reduced-state version of the original model by Jafri et al. (1998), which was implemented in the MATLAB code provided by the original author. The version of the model published in Saucerman et al. (2003) showed a complete version of the model as in Jafri et al. (1998), replicated to represent phosphorylated states.

3.2Parameter Values

The baseline parameter values are provided in the supplementary material of Saucerman et al. (2003) and in the MATLAB code provided by Saucerman. Missing values on cell geometry are provided in Table 1 which are used in all Figure reproductions.

NameValueSource
Myocyte volume36.8 pL36.8\ pLBers (2001)
Myocyte density1.06 g/mL1.06\ g/mLMendez et al. (1960)

Table 1:Parameters sourced for the model, which were not included in the primary paper.

It was found that the density of myocytes needed to increase by 6.4 times in order to better match the primary paper figures.

FigureParameter valuesNumber of iterationsT (min)
2a118
2b[Iso]=logspace(-3,0.45,8) μM[Iso^\dagger] = \texttt{logspace(-3,0.45,8)}\ \mu M6
[IBMX]=100 μM[IBMX] = 100\ \mu M8
2c, 2d[Iso]=1 μM[Iso] = 1\ \mu M61
[[PKI]={0,0.052} μM[[PKI] = \{0, 0.052\}\ \mu M
[cAMP]=logspace(-1.9,2,6) μM[cAMP^\dagger] = \texttt{logspace(-1.9,2,6)}\ \mu M
2e[Iso]=logspace(-6,0,20) μM[Iso^\dagger] = \texttt{logspace(-6,0,20)}\ \mu M2030
2f[Iso]={0,1} μM[Iso^\dagger] = \{0, 1\}\ \mu M21
2g, 2h[Iso]={0,1} μM[Iso^\dagger] = \{0, 1\}\ \mu M21
3a[Iso]=logspace(-3,0,8) μM[Iso^\dagger] = \texttt{logspace(-3,0,8)}\ \mu M327
[AC]tot=3[AC]_{tot} \asteq 3
[β1ARtot]=3[\beta1AR_{tot}] \asteq 3
KGsα/=3K_{Gs\alpha } \diveq 3
3b[Iso]=logspace(-3,0,8) μM[Iso^\dagger] = \texttt{logspace(-3,0,8)}\ \mu M327
[AC]tot =1.5[AC]_{tot}\ \asteq 1.5
[β1ARtot] =1.8[\beta1AR_{tot}]\ \asteq 1.8
KGsα /=1.9K_{Gs\alpha }\ \diveq 1.9
4a, 4b[IBMX]=250 μM[IBMX] = 250\ \mu M33
[Iso]={1,0,1} μM[Iso^\dagger] = \{1,0,1\}\ \mu M
khyd ={1,0.05,0.05}k_{hyd}^\dagger\ \asteq \{1,0.05,0.05\}
5[Iso]=logspace(-4,0,8) μM[Iso^\dagger] = \texttt{logspace(-4,0,8)}\ \mu M2420
[PKAI,tot] =0.05[PKA_{I,tot}]\ \asteq0.05
[PKItot] =0.05[PKI_{tot}]\ \asteq0.05
6a[Iso]={1,1,1} μM[Iso^\dagger] = \{1,1,1\}\ \mu M31
[PLBp]t=0 ={1,0,1}[PLBp^\dagger]^{t=0} \ \asteq \{1,0,1\}
[LCCap]t=0 ={1,0,1}[LCCap^\dagger]^{t=0} \ \asteq \{1,0,1\}
[LCCbp]t=0 ={1,0,1}[LCCbp^\dagger]^{t=0} \ \asteq \{1,0,1\}
6b[Iso]={0,1,0,1} μM[Iso^\dagger] = \{0,1,0,1\}\ \mu M41
fracPLB ={1,1,0,0}^\dagger\ \asteq \{1,1,0,0\}
fracPLBo ={1,1,0,0}^\dagger\ \asteq \{1,1,0,0\}

Table 2:Specific protocol followed in the reproduction of each figure, run for a span of T minutes. =,/=\asteq, \diveq indicate multiplication or division of parameter by the specified value respectively. \dagger denotes the variable value being changed with every iteration. One iteration constitutes running the entire model once. For multiple species being changed for one figure: for the nnth loop iteration, the value of the species is given by the nnth value in the sequence given in the sequence (Column 2) for all marked species. (3a, 3b): four separate iterations for every [Iso]: one control, one with changed [ACtot][AC_{tot}] only, one with changed [β1ARtot][\beta1AR_{tot}] only, and one with change KGsαK_{Gs\alpha } only. (5): three separate iterations for every [Iso]: one control, one with changed [PKAI,tot][PKA_{I,tot}] only, one with changed [PKItot][PKI_{tot}] only. patch clamping: Figures 2g-h, 6b follow Protocol 1; Figure 2f follows Protocol 2 (see Table 3).

See primary paper or source code for an explanation of all chemical abbreviations.

ProtocolDetails
1
Iapp={10 μAμF if mod(t+0.9,1)<=5e30 μAμF otherwiseI_{app} = \begin{cases} 10\ \frac{\mu A}{\mu F} & \ \text{if } \texttt{mod}(t+0.9, 1) <= 5e-3\\ 0\ \frac{\mu A}{\mu F} & \ \text{otherwise} \end{cases}
2
Vclamp={10 mV if 59.1<t<59.555 mV otherwise V_{clamp} = \begin{cases} -10\ mV & \ \text{if } 59.1<t<59.5\\ -55\ mV & \ \text{otherwise } \end{cases}
where Iapp=VclampVmRclampI_{app} = \frac{V_{clamp}-V_m}{R_{clamp}}

Table 3:Patch clamp protocols. Rclamp=0.02 ΩR_{clamp} = 0.02\ \Omega

Experimental conditions specific to each Figure are given in Table 2. For Figures 2f-h and 6b, patch clamp protocols are defined in Table 3.

This system has 44 dependent variables. Unless stated otherwise, all figure reproductions are subject to the same initial conditions (see source code). These values are embedded in the script main.m.

4Computational Simulation

The system of differential-algebraic equations is solved in MATLAB using stiff solver ode15s Shampine & Reichelt, 1997 with automatic/adaptive time stepping, RelTol=1e6RelTol = 1e-6, MaxStep=5e3MaxStep = 5e-3. No further packages are required; any machine with the full ODE suite (which includes ode15s solver) should be able to run all simulations.

Run the script main.m. The user will be prompted to specify the desired figure reproduction: input the appropriate string of format iaia where ii is the figure number and aa is the subfigure letter: for example, 2b2b.

Figures are highly recommended to be reproduced using parallel computing (parpool, >=4 workers) in MATLAB due to the size and stiffness of the code. By default, this is already set up in the script. Otherwise, computations will take up to 30 minutes per figure, on a 2.8 Ghz / 4 core machine running on Windows 10. Once running, keep parpool open (change parallel preferences: adjust time to automatically shut down pool after appropriate amount of time) such that time is not lost waiting for the pool to start up on every simulation.

FigureSpeciesOriginal valueValue modified for MATLAB
2bmax([Iso])10 μM10\ \mu M2.8 μM2.8\ \mu M
2c/d[PKI]60 nM60\ nM52 nM52\ nM
2fPatch clamp voltageholding 80 mV-80\ mV, step 0 mV0\ mVholding 55 mV-55\ mV, step 10 mV-10\ mV
4a[IBMX]1 mM1\ mM250 μM250\ \mu M
4bmax([Iso])10 μM10\ \mu M1 μM1\ \mu M
5PKAItot_{tot}05% of baseline
PKItot_{tot}05% of baseline
6a[Iso]-1 μM1\ \mu M
6b[Iso]-1 μM1\ \mu M
fracPLB-0

Table 4:Modifications made to stated values of primary paper, in the reproduction of the figures.

The original implementation of the model was performed using Berkeley-Madonna software, with a different ODE solver to MATLAB Saucerman et al., 2003. As some simulation values stated in the primary paper resulted in the model being unsolvable in the MATLAB routine, we had to make adjustments as outlined in Table 4.

5Model Results

We employ Engauge Digitizer software Mitchell et al., 2020 to extract data from Figures of Saucerman et al. (2003) to present alongside results of the present work.

The primary data (x) of Figure 2 with our reproduction of all subfigures (-). Note: (2c) PKA\ activation = \frac{[PKACI]}{2*[PKAI_{tot}]}

Figure 1:The primary data (xx) of Figure 2 with our reproduction of all subfigures (-). Note: (2c) PKA activation=[PKACI]2[PKAItot]PKA\ activation = \frac{[PKACI]}{2*[PKAI_{tot}]}

The primary data (x) of Figure 3 with our reproduction of all subfigures (-). Note: lines in (3b) showing modifications to [AC]_{tot}, \beta 1AR, K_{Gs\alpha } have been grouped into an average ‘augmentation’ line.

Figure 2:The primary data (xx) of Figure 3 with our reproduction of all subfigures (-). Note: lines in (3b) showing modifications to [AC]tot,β1AR,KGsα[AC]_{tot}, \beta 1AR, K_{Gs\alpha } have been grouped into an average ‘augmentation’ line.

The primary data (x) of Figure 4 with our reproduction of all subfigures (-). (4a): primary data shown in blue, reproduction in red. AC\ activity=cAMP_{tot}. (4b): f_{cAMP}=\frac{cAMP_{tot}}{cAMP_{tot,max}}. Annotated numbers on the graph refer to [Iso] set for that experiment.

Figure 3:The primary data (xx) of Figure 4 with our reproduction of all subfigures (-). (4a): primary data shown in blue, reproduction in red. AC activity=cAMPtotAC\ activity=cAMP_{tot}. (4b): fcAMP=cAMPtotcAMPtot,maxf_{cAMP}=\frac{cAMP_{tot}}{cAMP_{tot,max}}. Annotated numbers on the graph refer to [Iso] set for that experiment.

The primary data (x) of Figure 5 with our reproduction (-).

Figure 4:The primary data (xx) of Figure 5 with our reproduction (-).

(6a) The primary data (x) of Figure 6a with our reproduction of all subfigures (-). (6b) Primary data shown in blue, reproduction in red.

Figure 5:(6a) The primary data (xx) of Figure 6a with our reproduction of all subfigures (-). (6b) Primary data shown in blue, reproduction in red.

The reproduction of all figures of Saucerman et al. (2003) are given in Figures 1-5, which show good agreement against the data of the primary paper. Solid lines show the output from this curation effort, and crosses show discrete points found by the Engauge Digitizer of the figures originally published.

6Discussion

We demonstrate that figures of Saucerman et al. (2003) can be reproduced given some additions (Table 1) and adjustments (Table 4) made to model parameter values. There are slight deviations of the curated model output from the original data which can be attributed to parameter modifications we had to make to ensure a simulation runs to completion (see Figures 1c, 1d, 3a, 3b). Other mismatches may be owing to differences from the original implementation in Berkeley-Madonna, which had different solver settings to this effort. The original code is no longer available.

As this paper follows the tenets of FAIR Wilkinson et al., 2016, further credibility is lent to the original model. The model has already proven to be valuable to the mathematical modelling community in understanding the workings of β\beta-adrenergic receptor signalling within cardiac myocytes.

References
  1. Saucerman, J. J., Brunton, L. L., Michailova, A. P., & McCulloch, A. D. (2003). Modeling β-adrenergic control of cardiac myocyte contractility in silico. Journal of Biological Chemistry, 278(48), 47997–48003.
  2. MATLAB. (2017). version 9.3.0713579 (R2017b). The MathWorks Inc.
  3. Cuellar, A. A., Lloyd, C. M., Nielsen, P. F., Bullivant, D. P., Nickerson, D. P., & Hunter, P. J. (2003). An Overview of CellML 1.1, a Biological Model Description Language. SIMULATION, 79(12), 740–747. 10.1177/0037549703040939
  4. Jafri, M. S., Rice, J. J., & Winslow, R. L. (1998). Cardiac Ca2+ dynamics: the roles of ryanodine receptor adaptation and sarcoplasmic reticulum load. Biophysical Journal, 74(3), 1149–1168.
  5. Bers, D. (2001). Excitation-contraction coupling and cardiac contractile force (Vol. 237). Springer Science & Business Media.
Funding and support by:
International Union of Physiological SciencesAuckland Bioengineering InstituteDigital ScienceVPH Institute