# A kinetic model of β-adrenergic control in cardiac myocytes

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.

## 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 ($\beta$1AR), 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:

$\beta$1AR 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.

Name | Value | Source |
---|---|---|

Myocyte volume | $36.8\ pL$ | Bers (2001) |

Myocyte density | $1.06\ g/mL$ | Mendez 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.

Figure | Parameter values | Number of iterations | T (min) |
---|---|---|---|

2a | 1 | 18 | |

2b | $[Iso^\dagger] = \texttt{logspace(-3,0.45,8)}\ \mu M$ | 6 | |

$[IBMX] = 100\ \mu M$ | 8 | ||

2c, 2d | $[Iso] = 1\ \mu M$ | 6 | 1 |

$[[PKI] = \{0, 0.052\}\ \mu M$ | |||

$[cAMP^\dagger] = \texttt{logspace(-1.9,2,6)}\ \mu M$ | |||

2e | $[Iso^\dagger] = \texttt{logspace(-6,0,20)}\ \mu M$ | 20 | 30 |

2f | $[Iso^\dagger] = \{0, 1\}\ \mu M$ | 2 | 1 |

2g, 2h | $[Iso^\dagger] = \{0, 1\}\ \mu M$ | 2 | 1 |

3a | $[Iso^\dagger] = \texttt{logspace(-3,0,8)}\ \mu M$ | 32 | 7 |

$[AC]_{tot} \asteq 3$ | |||

$[\beta1AR_{tot}] \asteq 3$ | |||

$K_{Gs\alpha } \diveq 3$ | |||

3b | $[Iso^\dagger] = \texttt{logspace(-3,0,8)}\ \mu M$ | 32 | 7 |

$[AC]_{tot}\ \asteq 1.5$ | |||

$[\beta1AR_{tot}]\ \asteq 1.8$ | |||

$K_{Gs\alpha }\ \diveq 1.9$ | |||

4a, 4b | $[IBMX] = 250\ \mu M$ | 3 | 3 |

$[Iso^\dagger] = \{1,0,1\}\ \mu M$ | |||

$k_{hyd}^\dagger\ \asteq \{1,0.05,0.05\}$ | |||

5 | $[Iso^\dagger] = \texttt{logspace(-4,0,8)}\ \mu M$ | 24 | 20 |

$[PKA_{I,tot}]\ \asteq0.05$ | |||

$[PKI_{tot}]\ \asteq0.05$ | |||

6a | $[Iso^\dagger] = \{1,1,1\}\ \mu M$ | 3 | 1 |

$[PLBp^\dagger]^{t=0} \ \asteq \{1,0,1\}$ | |||

$[LCCap^\dagger]^{t=0} \ \asteq \{1,0,1\}$ | |||

$[LCCbp^\dagger]^{t=0} \ \asteq \{1,0,1\}$ | |||

6b | $[Iso^\dagger] = \{0,1,0,1\}\ \mu M$ | 4 | 1 |

fracPLB$^\dagger\ \asteq \{1,1,0,0\}$ | |||

fracPLBo$^\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 $n$th loop iteration, the value of the species is given by the $n$th 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 $[AC_{tot}]$ only, one with changed $[\beta1AR_{tot}]$ only, and one with change $K_{Gs\alpha }$ only. (5): three separate iterations for every [Iso]: one control, one with changed $[PKA_{I,tot}]$ only, one with changed $[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.

Protocol | Details |
---|---|

1 | $I_{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 | $V_{clamp} =
\begin{cases}
-10\ mV & \ \text{if } 59.1<t<59.5\\
-55\ mV & \ \text{otherwise }
\end{cases}$ |

where $I_{app} = \frac{V_{clamp}-V_m}{R_{clamp}}$ |

Table 3:Patch clamp protocols. $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 = 1e-6$, $MaxStep = 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 $ia$ where $i$ is the figure number and $a$ is the subfigure letter: for example, $2b$.

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.

Figure | Species | Original value | Value modified for MATLAB |
---|---|---|---|

2b | max([Iso]) | $10\ \mu M$ | $2.8\ \mu M$ |

2c/d | [PKI] | $60\ nM$ | $52\ nM$ |

2f | Patch clamp voltage | holding $-80\ mV$, step $0\ mV$ | holding $-55\ mV$, step $-10\ mV$ |

4a | [IBMX] | $1\ mM$ | $250\ \mu M$ |

4b | max([Iso]) | $10\ \mu M$ | $1\ \mu M$ |

5 | PKAI$_{tot}$ | 0 | 5% of baseline |

PKI$_{tot}$ | 0 | 5% of baseline | |

6a | [Iso] | - | $1\ \mu M$ |

6b | [Iso] | - | $1\ \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 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.

- 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. - MATLAB. (2017).
*version 9.3.0713579 (R2017b)*. The MathWorks Inc. - 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 - 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. - Bers, D. (2001).
*Excitation-contraction coupling and cardiac contractile force*(Vol. 237). Springer Science & Business Media.