Mathematical modeling improves EC50 estimations from classical dose–response curves

(1)

Mathematical modeling improves EC50

estimations from classical dose–response curves

Elin Nyman, Isa Lindgren, William Lövfors, Karin Lundegård, Ida Cervin, Theresia Arbring, Jordi Altimitas and Gunnar Cedersund

Linköping University Post Print

N.B.: When citing this work, cite the original article.

Original Publication:

Elin Nyman, Isa Lindgren, William Lövfors, Karin Lundegård, Ida Cervin, Theresia Arbring, Jordi Altimitas and Gunnar Cedersund, Mathematical modeling improves EC50 estimations from classical dose–response curves, 2015, The FEBS Journal, (282), 5, 951-962.

http://dx.doi.org/10.1111/febs.13194 Copyright: Wiley: FEBS Journal

http://onlinelibrary.wiley.com/journal/10.1111/%28ISSN%291742-4658

Postprint available at: Linköping University Electronic Press

(2)

Mathematical modeling improves EC

50

estimations from classical

dose-response curves

Elin Nyman1,2, Isa Lingren3, William Lövfors1, Karin Lundengård1,4, Ida Cervin1, Theresia Arbring Sjöström1,5, Jordi Altimiras3, Gunnar Cedersund1,6

1

Integrative Systems Biology, Department of Biomedical Engineering, Linköping University, SE-581 85 Linköping, Sweden

2

CVMD iMED DMPK AstraZeneca R&D, SE-431 83 Mölndal, Sweden 3

AVIAN Behavioural Genomics and Physiology, Department of Physics, Chemistry and Biology, Linköping University, SE-581 83 Linköping, Sweden

4

Radiological Sciences, Department of Medicine and Health, Linköping University, SE-581 85 Linköping, Sweden

5

Laboratory of Organic Electronics, Department of Science and Technology, Linköping University, SE-601 74 Norrköping, Sweden

6

Cell biology, Department of Clinical and Experimental Medicine, Linköping University, SE-581 85 Linköping, Sweden

Corresponding author:

Gunnar Cedersund, Department of Clinical and Experimental Medicine, Linköping University, SE58185 Linköping, Sweden, Email: gunnar.cedersund@liu.se, Phone: +46-702-512323, Fax: +46-10 1034149

(3)

The beta-adrenergic response is impaired in failing hearts. When studying beta-adrenergic function in vitro, the half-maximal effective concentration (EC50) is an important measure of ligand response. We previously measured the in vitro contraction force response of chicken heart tissue to increasing concentrations of adrenaline, but noted a decreasing response at high

concentrations. The classical interpretation of such data would be to assume maximal response before the decrease and fit a sigmoid curve to the remaining data to determine EC50. We have applied a mathematical modeling approach to instead interpret the full dose-response curve in a new way. The developed model predicts a non-steady-state caused by short resting time between increased concentrations of agonist, which affect the dose-response characterization. Therefore, an improved estimate of EC50 can be calculated using steady-state simulations of the model. The model-based estimation of EC50 is further refined with additional time-resolved data to decrease the uncertainty of the prediction. The resulting model-based EC50 (180-525 nM) is higher than the classical interpreted EC50 (46-191 nM). Mathematical modeling thus makes it possible to reinterpret already obtained datasets, and to make accurate estimates of EC50 even when steady-state measurements are not experimentally feasible.

adrenaline / cardiac beta-adrenergic signaling / dynamic mathematical modeling / EC50 / ordinary differential equations

(4)

Introduction

Normal cardiac function is heavily dependent on adrenergic and noradrenergic stimulation of beta-adrenergic receptors (βARs) in the myocardium. In diseased hearts, the βAR function is often impaired and βARs are therefore a classic and well-studied target for medical intervention. In basic research, in vitro cardiac βAR function is commonly measured as the contractile

response of heart tissue to βAR stimulation by bioamines such as adrenaline (ADR) and noradrenaline (binds to α and βARs) or the structurally similar synthetic agonist isoproterenol (ISO, βAR-specific).

The βAR pathway is initiated by the activation of a stimulatory G-protein that in turn activate adenylyl cyclase (AC) to produce the second messenger cyclic AMP (cAMP) (reviewed in Rockman, et al. [1]). cAMP has many downstream targets, but one of the main ones is protein kinase A (PKA), which in turn phosphorylates multiple targets, including the βAR itself leading to βAR uncoupling from its G-protein, and thereby to desensitization [2]. The βAR pathway is regulated through several layers of feedback mechanisms for fine-tuning and protection from over-stimulation. Like PKA, other enzymes such as phosphodiesterase (PDE) and G-protein receptor kinases (GRK), can through different mechanisms negatively affect the βAR response [3]. Desensitization of the βAR system, possibly also leading to full downregulation by

sequestration and degradation of the receptor, is crucial in protecting the heart from long-term adrenergic over-stimulation [3]. Dysregulation of the βAR system can be the foundation for, or the result of, different cardiovascular diseases and is of great medical importance [4-8].

Dose-response curves and the calculated half maximal effective concentration (EC50) are often used as a method to measure drug effects in biological systems. The calculation of EC50 rely on a number of assumptions such as i) a sigmoidal response to increasing agonist concentrations, and ii) a reached response “plateau” (steady-state) before the next cumulative drug concentration is added to the system. Due to experimental limitations, these criteria may not always be fulfilled. To complicate matters further, scientific reports sometimes only present the calculated EC50 and not the raw data for the dose-response curve. This limits the conclusions one can draw from the

(5)

analyzing dose-response curves may appear straight-forward, but can also be complicated to the point that information is masked and the interpretation of EC50 is confounded.

We have previously studied the βAR response to ADR in chicken ventricular tissue [9] and observed a decreased contraction response at higher agonist concentrations (unpublished data). This decrease has been related to receptor desensitization/downregulation [10-12] and has also been ascribed to the activation of inhibitory G-protein dependent pathways [2, 13, 14]. In this study, we use mathematical modeling to reinterpret adrenergic dose-response curves that display a decrease in the response to high concentrations of agonist. Through this modeling approach, we show an alternative, new explanation to the drop in response to high concentrations of agonist as a consequence of non-steady-state measurements, which potentially confounds the interpretation of previously calculated EC50 values. Because of the time-limited explant viability, the

experimental protocol for generating the initial dose-response curves did not allow longer intervals between agonist doses to obtain proper steady-state measurements. As an alternative, we suggest that such experimental limitations can be circumvented through steady-state model simulations based on the initial data.

(6)

Results

The dose-response curve exhibits a decrease in contraction for high concentrations of adrenaline

We stimulated chicken ventricular tissue with increasing concentrations of ADR according to an experimental protocol (Figure 1A) to determine the EC50. The force of contraction in the cardiac tissue samples increased with increased concentrations of ADR, up to a maximal value of 50 % increase over basal force of contraction (Figure 1B). The increase in the force of contraction continued until 1000 nM of ADR was added, but with 3000 nM the increase changed to a decrease in force of contraction (Figure 1B). This behavior was clearly seen in 4 out of 6 individual experiments (Figure 1C, black), and in the 2 remaining experiments there are more variation in the response (Figure 1C, yellow and orange). If one assumes that maximal activation is obtained at 1000 nM and fit a sigmoidal curve to the data (Figure 1D), the calculated

confidence interval of EC50 is 46-191 nM (Figure 1E) (Materials and Methods). However, we sought for a sounder way of interpreting these data using mathematical modeling.

A minimal mathematical model gives new insights

The obtained data (Figure 1B) contains dynamic information which is best utilized with dynamic mathematical modeling frameworks. We therefore chose to develop a dynamic model based on ordinary differential equations (ODEs). In the design of the ODE model, we considered both our obtained data and available knowledge about βAR signaling. βAR signaling is complex, with many proteins involved and several layers of feedbacks, and to fully account for this complexity in the model, we would need additional data for the intracellular intermediaries. Instead, we developed a minimal model for βAR receptor signaling, i.e. a model with as few states as possible while still agreeing with available data. The developed minimal model (Figure 2A) contains 4 states, including cyclic AMP (cAMP) and the βAR receptor in three states: basal inactive (β), active (βact), and phosphorylated inactive (βphos), 7 parameters

(kbasal,kH,k1,k2,k3,k4,cAMP0), and one input signal (H). The input signal

corresponds to the scheme of added ADR and ISO concentrations in the experimental setup (cf. Figure 1A). With these notations, the corresponding ODEs for the minimal model are:

(7)

d/dt(βact) = (kbasal+kH·H)·β – k2·βact d/dt(βphos) = k2·βact - k1·βphos

d/dt(cAMP) = cAMP0 + k3·βact – k4·cAMP

The model thus allow for a basal activation of the receptors (kbasal·β) with an increase in activation with addition of the agonist (kH·H·β). The activated receptor is phosphorylated before returning to the inactive state where new activation can occur. The amount of cAMP is changed both with a constant positive term (cAMP0), with a term that is proportional to the amount of active receptors (k3·βact), and with a term that gives a breakdown proportional to the amount of cAMP (k4·cAMP). As model output, we use the increase in cAMP scaled with a parameter (kscale) and compare with the experimentally measured force of contraction.

y = kscale·(cAMP – cAMP(0))

The scaling parameter (kscale) was set to 1 to reduce the number of parameters to estimate. Initial conditions were not known; we therefore used 100 % as the total number of receptors, and an arbitrary chosen value of 10 for cAMP:

β(0) = 100 βact(0) = 0 βphos(0) = 0 cAMP(0) = 10

The model was simulated for 1000 minutes without input (H = 0) to obtain steady-state values for the model states for each new set of parameters that was tested. These parameter-dependent steady-state values were used as initial conditions in the simulated experiments. In the simulated experiments, the value of the input was changed according to the experimental protocol (cf. Figure 1A).

The parameters of the model were estimated to fit with the data in Figure 1B. The best fit model simulations and the experimental dose-response curve are in good agreement (Figure 2B, bold)

(8)

(χ2

(minimal model)=60.5 < 65). We searched not only for the best fit parameter sets, but also for all parameter sets that provide model simulations with good enough agreement with

experimental data according to χ2-statistics (Materials and Methods). An approximation of all found acceptable parameters demonstrate the different possible behaviors with statistically acceptable fit (Figure 2B, lines). These simulations all show the decrease in the dose-response curve at high concentrations of ADR (Figure 2B, lines). Time-resolved simulations with continuous responses for all doses of ADR demonstrate dynamic overshoots as the mechanism behind the decrease (Figure 3A). The predicted dynamic overshoots give accumulating responses during the simulation and provide the agreement with data (Figure 3B). This prediction means that the decrease in contractility at the highest ADR concentrations is a dynamic phenomenon caused by not waiting long enough until the next dose is applied. The dynamic overshoot is in the model created by a slow return of the phosphorylated inactive receptors to the basal state where new activation can occur. Using the acceptable parameters in Figure 2B, the dynamic overshoot behavior is predicted for the whole range of parameters for a high concentration of ADR (1000 nM) (Figure 3C).

Data contains uncertainties that aggravate a precise prediction of EC50

With the developed and parameterized model it is possible to perform simulations of a virtual steady-state experiment. Using the acceptable parameters of the model, we create steady-state stimulations (Figure 3D) and corresponding dose-response curves (Figure 3E). However, the uncertainty of the data (reflected by the acceptable parameters) is large and the EC50 can only be predicted with a corresponding large uncertainty (44-369 nM).

Independent time-resolved data confirms the overshoot behavior

To obtain a better estimate of the EC50 we gathered new experimental data using a different experimental protocol. We performed a time-resolved experiment with the alternative,

structurally similar, agonist ISO. To maximize the information in the data and thus obtain a more precise parameterization of the model, we stimulated the cardiac tissue twice with an in between washing step (Figure 4A). The resulting force of contraction displayed a dynamic overshoot behavior and a slight desensitization with a lower response to the second stimuli (Figure 4A).

(9)

This result confirmed the overshoot prediction from the minimal model and thus strengthens the minimal model as a valid description of the system.

The same minimal model can explain the time-resolved data

Using the same developed minimal model (Figure 2A), we searched for parameters with an acceptable agreement to the time-resolved, ISO stimulated data, but found no such parameters (χ2(minimal model)=44.6 > 31). The reason for the bad agreement was mainly the first data point after addition of ISO (Figure 4B, red circle), where the response to the stimulation is slower than the model can account for (Figure 4B, red line). The underlying reason for the slow response may be that it takes time for the added ISO to diffuse in the solution and the tissue. This issue arises since we increase the concentration of the agonist directly from 0 to 100 nM instead of the step-wise increase in the dose-response measurement. We handled this problem with an extra time-delay parameter that gives freedom to the model for a later effect of the added ISO, and found that such a parameter substantially increased the agreement with data (Figure 4B, green line) (χ2(minimal model)=18.7 < 31). We performed a profile likelihood analysis [15] for the time-delay parameter (Figure 4C), and found that 0.65 minutes (39 seconds) is the optimal value of the time-delay and we chose that value for further analysis.

Accounting for the different agonists using the minimal model

To obtain a better prediction of EC50, we wanted to simultaneously use the datasets for ADR stimulated dose-response and ISO stimulated time-response. To do so, we introduced a single parameter for the ratio of the difference in affinity between ADR and ISO (kratio). This parameter was set to 1 when fitting to the ADR data and allowed to change when fitting to the ISO data. In other words, the new equations are:

d/dt(β) = -(kbasal+kratio·kH·H)·β + k1·βphos d/dt(βact) = (kbasal+kratio·kH·H)β – k2·βact

d/dt(βphos) = k2·βact - k1·βphos

(10)

We performed a profile likelihood analysis [15] for the parameter kratio and the 95 % confidence interval for the parameter is found for ratios between 3 to 9 (Figure 5). This value represent that ISO has 3-9 fold higher affinity for the βAR receptor than ADR, which is in agreement with the literature [16-18]. This agreement again serves as a confirmation that the developed minimal model is a valid description of the studied system. We thus chose

3<kratio<9 as boundaries in the following analyses.

A model-based estimation of EC50

The combined dataset (ADR stimulated dose-response and ISO stimulated time-response) should contain more information than the original dataset, and therefore allow for a more well-defined estimation of the EC50 for ADR stimulated force of contraction in cardiac tissue. We gathered parameters with acceptable fit to the combined dataset (Table I, Figure 6A). Using these parameters, we simulated an experiment where steady-state is obtained between every increase in ADR concentration (Figure 6B). The corresponding dose-response curve for the steady-state responses gives EC50 of 193-493 nM ADR for all acceptable parameter sets (Figure 6C). We also searched specifically for parameters that stretch this prediction, that is parameters that give the minimal and maximal EC50 while still agreeing with all data (see Materials and Methods) and the possible range of parameters was extended to a calculated 95 % confidence interval of 180-525 nM with this extended search (Figure 6D). The difference between the straight-forward

calculated confidence interval of EC50 (46-191 nM), and the EC50 obtained with the model-based approach (180-525 nM) is substantial (Figure 6E). However, the confidence intervals are not completely separated, but have a slight overlap of <10 % of the interval lengths.

(11)

Discussion

The calculation of EC50 is one of the most common calculations in biological and medical research. We show herein that this calculation is not always straight-forward and that it is important to present not only the calculated EC50 value, but also the data behind the calculation, that is the dose-response curve. We present a mathematical modeling approach for analysis of EC50 in dose-response experiments in cardiac tissue. With this approach, we provide evidence for a novel interpretation for a decrease of contraction at high concentrations of ADR. This interpretation implies that steady-state was not reached in the experimental setup with the decrease in contractions at high concentrations of ADR as a consequence. We also show how to use the model-based interpretation to recalculate the EC50 for ADR stimulation of contraction. Using our approach, it is thus possible to “rescue” already obtained experimental data for this system.

The apparent solution to the non-steady-state measurements performed herein would be to redo the experiments with longer resting times between increased doses of ADR in the experimental setup. However, this is not possible since the contractile performance of heart tissue samples deteriorates over time. Instead, we employ mathematical modeling to simulate the wanted

experimental conditions, where steady state in the model is obtained before new doses of ADR is added. In fact, mathematical modeling is necessary to obtain a sound interpretation of data for this experimental setup, where we have sensitive tissue samples that do not allow for steady-state measurements. The developed mathematical model can also save experimental time and work in future experiments.

The main strength of the mathematical modeling method used herein is the use of minimal models in combination with evaluation of all found parameter sets that are in agreement with available data [19-21]. This is important in cases where model parameters cannot be uniquely estimated due to uncertainty and sparseness of data. Using all acceptable parameter sets, general model behavior can be studied that do not rely on the specific values of the parameters. For further analysis of the properties of the acceptable parameters, their co-variation structure can be evaluated with principal component analysis [22]. Such evaluation increases the understanding

(12)

of model properties and can be used for model reduction. Another similar analysis that is useful in cases of parameter unidentifiability combines the Fisher Information Matrix with profile likelihood analysis to semi-automatically decide identifiable parameter combinations [23].

Mathematical modeling has previously been used to understand different aspects of βAR receptor signaling. Saucerman, et al. [24] developed a comprehensive mathematical model that consists of the βAR receptor signaling chain from receptor activation to cell contraction. This model has been further extended in several steps [25-27] and analyzed and reduced down to a minimal model with 4 states [28]. Bondarenko [29] recently presented a compartmentalized mathematical model of βAR signaling. Other models of interest concern an increased

understanding of the mechanisms behind βAR desensitization, since desensitization give rise to overshoots in data. Xin, et al. [30] have evaluated the different contributions of GRK and PDE signaling to overshoot dynamics in cAMP data and found that inhibition of PDE as well as combined PKA and GRK inhibition removed the overshoot in a cell line. Violin, et al. [3] studied the same cell system and found that direct receptor inactivation by PKA is non-important, but that GRK plays an important role for the overshoot behavior. Shankaran, et al. [31] performed a theoretical study of G-protein coupled receptor desensitization and compared with receptor downregulation through internalization and found that both mechanisms provide efficient input-output coupling. None of the mentioned mathematical models were developed to understand the underlying mechanisms to a decrease in the dose-response curve during

stimulation.

The observed overshoot behavior is found not only in βAR signaling, but also in other signaling networks such as the insulin signaling network. In the insulin signaling network, both the insulin receptor and the insulin receptor substrate respond to insulin with a dynamic overshoot [20, 32]. This insulin signaling overshoot seems to spread to some, but not all, branches of the signaling network [33]. The dynamic overshoot behavior in the insulin signaling network arises from negative feedback mechanisms and act to decrease the cellular response to a stimulus for high concentrations [20].

(13)

In this study, we use the change in cAMP concentration in the mathematical model as a direct readout for the experimentally measured force of contraction. The force of contraction has been shown to correlate with cAMP levels during ADR stimulation in perfused rat hearts [34]. Nevertheless, this readout is an over-simplification. The main reason for the use of this

simplification is that we aim to obtain a minimal model of the system to be able to analyze the mechanisms behind the decrease in the dose-response curve in a conclusive manner. The model is also a simplification when it comes to the βAR signaling pathway. To keep the model

minimal, we have not included the difference between β1 and β2 adrenergic receptors with the resulting cross-signaling via inhibitory G-proteins [14, 35]. Finally, we have simplified the model when it comes to negative feedback mechanisms. In βAR receptor signaling, there are several levels of negative feedback mechanisms: PKA- and GRK-mediated inhibitory

phosphorylation of receptors [2, 14, 26, 30], and internalization of receptors to the cell interior where activation cannot occur [36]. We incorporated a single negative feedback mechanism: inactivation due to inhibitory phosphorylation of the receptors, in the developed mathematical model. However, we expect the conclusions herein to also hold for a more complex/realistic mathematical model of βAR signaling.

The use of different agonists for the dose-response and the time-resolved experiments may not be ideal; ADR and ISO are structurally similar and both stimulate βAR receptors, but have different βAR affinity and efficacy. However, we are able to demonstrate the possibility to combine data from pre-existing datasets where different agonists have been used and thereby demonstrate the strength of mathematical modeling in a G-protein coupled receptor system. Furthermore, the model-based estimation of the difference in affinity between ADR and ISO (3-9 times higher affinity for ISO) is in agreement with literature [16-18] and this result confirms that the developed minimal model is detailed enough to capture the response to both ADR and ISO stimulation.

In the time-resolved data, we obtained a rebound effect in the washing step, namely that the removal of the stimulus gives a lower than basal response of the system (Figure 4A). Such a rebound effect is unwanted in drug dosing and evaluation of treatment [37]. The rebound effect is also explained with the same overshoot mechanism in the mathematical model. This means

(14)

that the developed model can potentially be used to study rebound effects and the mechanisms behind it, and importantly, to study how to remove such unwanted effects.

In summary, we have shown how mathematical modeling sometimes is necessary to calculate as simple and common biological measures as EC50, and exemplified the model-based approach on the well-studied and medically important βAR signaling system in cardiac tissue.

(15)

Materials and Methods

Experimental methods

The methods have been described elsewhere [9]. Briefly, trabeculae from 5 week old chickens were isolated and suspended between an immobile rod in one end and a force transducer in the other. The tissue was submerged in 37°C modified Ringer buffer solution (composition in mM: 138 NaCl, 3 KCl, 3 CaCl2, 1.8 MgCl2, and 10 HEPES, 5 Sodium Pyruvate, Tris to pH 7.4) bubbled with 100% O2. Platinum electrodes on either side of the tissue stimulated continuous contractions at a frequency of 1.5 Hz, pulse duration 20 ms and a voltage of 10-20 V. After ~60 min of stable contractions, the muscle tension was set to that which produced 80% of the

maximal twitch force according to the Frank-Starling law and was kept at this tension throughout the experiment. A 5 minutes sequence of baseline contractions was recorded before any drugs were added. In the dose-response experiment, doses of ADR (0.01 – 10000 nM) were added in a cumulative fashion every 150 second. The contraction force was measured as the average of 20 contraction cycles within the last 20 seconds before the addition of the next concentration of agonist. In the time-resolved experiment, each trabecula was stimulated with 100 nM ISO for 10 minutes. The stimulation was interrupted by washing the trabecula and allowing it to rest for 5 minutes before another 10 minute stimulation with 100 nM ISO. The mean force of contraction was calculated for each minute of the recorded data.

Classical EC50 estimation

We fitted sigmoidal standard curves to the first 7 data points in Figure 1D. We used the 4 parameter logistic equation:

𝐹(𝑐𝑜𝑛𝑐) = 𝐴 + (𝐷 – 𝐴) / (1 + 10^(log (𝐶) − log (𝑐𝑜𝑛𝑐))·𝐵) where conc is the concentration, A is the asymptote value for low concentrations, D is the asymptote value for high concentrations, B is the steepness of the curve, and C is the inflection point. We chose best fit values A = 0.983, D = 1.48, and B = 2, and varied the values of C and calculated the corresponding fit for the sigmoidal standard curve obtained using χ2 statistics (see

(16)

below). Regarding degrees of freedom, we used n=7 data points minus 4 parameters, i.e. χ2 (7-4,0.05)=7.8.

Mathematical modeling

The developed minimal mathematical models are based on ordinary differential equations (ODEs):

𝑥̇ = 𝑓(𝑥, 𝑝, 𝑢) 𝑦 = 𝑔(𝑥, 𝑝, 𝑢)

where x represent the states of the model, p represents the parameters of the model, u represents the model inputs, y represent the output of the model that correspond to experimentally measured responses, and f and g are nonlinear functions that describe the mechanistic assumptions of the model. The ODEs were simulated using the Systems Biology Toolbox for Matlab [38, 39]. The formulation of ODEs from an interaction-graph is described in [40] and in the supplementary material of [21]. Mass-action kinetics was used in the rate-equations to obtain a minimal number of parameters. These parameters were estimated using our modified version of the function simannelingSBAO in Systems Biology Toolbox. The modified version is complemented with multiple downhill simplexes for a more comprehensive search of the parameter space [19, 20]. All Matlab model files and scripts required to simulate the figures herein are available in the catalogue “ModelFiles.zip”.

Statistical analyses

To evaluate the agreement between model simulations and experimental data we calculated the

cost, i.e. the sum of squares of the residuals weighted by the variance of the data:

𝑐𝑜𝑠𝑡(𝑝) = ∑ (𝑦𝑖 − 𝑦̂(𝑝)𝑖 𝑆𝐸𝑀_𝑖 )

2 𝑁

(17)

where N is the number of data points, yi – ŷ(p)i is the ith residual between the data point and the

simulated model response normalized by the variance (SEM). The simulated model response

ŷ(p)i depends on the values of the parameters (p). The calculated sum follows a χ2 distribution.

Therefore, χ2-tests were used to decide upon whether models should be rejected [41]. In χ2-tests, the level of significance and the degrees of freedom decide the border of rejection from a χ2 distribution. We used a significance level of 0.05 and degrees of freedom where taken as the number of data points minus 1 parameter for normalization minus the number of parameters in the model (7, 8, or 9). For the dose-response curve with n=56 data points the border of rejection is χ2

(55-1-7,0.05)=64, for the time-response curve with n=29 data points the border of rejection is χ2

(29-1-8,0.05)=31, and for the combined data set with n=56+29 data points the border of rejection is χ2

(85-2-9,0.05)=95. We saved all parameters that passed the χ2-test and picked out the parameter sets containing a maximal or minimal value of any of the parameters (Table I, column 3 and 4) to use as an approximation of the total set. In the approximation of the total set, we also picked the parameters with the best and worst agreement with data from the saved parameters. We used this approximation of the acceptable parameters in the figures. We performed profile likelihood calculations [15] for a few of the model parameters that were of extra interest for the analyses. We also performed parameter searches in the direction of specific predictions according to the methods developed in [19, 42]. In these searches, the parameters are forced in the direction of the prediction of interest, and the agreement with data is a side

condition that also should be fulfilled. These parameter searches are valuable for complex estimation problems where it is hard to find a good approximation for the total set of acceptable parameters.

(18)

Acknowledgements

We thank students from the course TSRT17 and TBMT19 for contributing in this project. This work was supported by LIST, Östergötland County Council, and the Swedish Research Council.

Author Contribution

JA and GC designed the project. EN and GC designed the mathematical modeling. EN, WL, KL, IC, and TS performed the mathematical modeling. IL and JA designed the experimental

procedures. IL and KL performed the experimental procedures. EN, IL, JA, and GC wrote the manuscript.

Conflict of Interest

The authors declare no conflict of interest.

(19)

References

1. Rockman, H. A., Koch, W. J. & Lefkowitz, R. J. (2002) Seven-transmembrane-spanning receptors and heart function, Nature. 415, 206-12.

2. Martin, N. P., Whalen, E. J., Zamah, M. A., Pierce, K. L. & Lefkowitz, R. J. (2004) PKA-mediated phosphorylation of the beta1-adrenergic receptor promotes Gs/Gi switching, Cell Signal. 16, 1397-403. 3. Violin, J. D., DiPilato, L. M., Yildirim, N., Elston, T. C., Zhang, J. & Lefkowitz, R. J. (2008) beta2-adrenergic receptor signaling and desensitization elucidated by quantitative modeling of real time cAMP dynamics, J Biol Chem. 283, 2949-61.

4. El-Armouche, A. & Eschenhagen, T. (2009) Beta-adrenergic stimulation and myocardial function in the failing heart, Heart Fail Rev. 14, 225-41.

5. Eschenhagen, T. (2008) Beta-adrenergic signaling in heart failure-adapt or die, Nat Med. 14, 485-7. 6. Sucharov, C. C. (2007) Beta-adrenergic pathways in human heart failure, Expert Rev Cardiovasc Ther.

5, 119-24.

7. Nikolaev, V. O., Moshkov, A., Lyon, A. R., Miragoli, M., Novak, P., Paur, H., Lohse, M. J., Korchev, Y. E., Harding, S. E. & Gorelik, J. (2010) Beta2-adrenergic receptor redistribution in heart failure

changes cAMP compartmentation, Science. 327, 1653-7.

8. Liggett, S. B., Cresci, S., Kelly, R. J., Syed, F. M., Matkovich, S. J., Hahn, H. S., Diwan, A., Martini, J. S., Sparks, L., Parekh, R. R., Spertus, J. A., Koch, W. J., Kardia, S. L. & Dorn, G. W., 2nd (2008) A GRK5 polymorphism that inhibits beta-adrenergic receptor signaling is protective in heart failure, Nat

Med. 14, 510-7.

9. Lindgren, I. & Altimiras, J. (2009) Chronic prenatal hypoxia sensitizes beta-adrenoceptors in the embryonic heart but causes postnatal desensitization, Am J Physiol Regul Integr Comp Physiol. 297, R258-64.

10. Lindgren, I. & Altimiras, J. (2013) Prenatal hypoxia programs changes in beta-adrenergic signaling and postnatal cardiac contractile dysfunction, Am J Physiol Regul Integr Comp Physiol. 305, R1093-101. 11. Rohlicek, C. V., Viau, S., Trieu, P. & Hebert, T. E. (2005) Effects of neonatal hypoxia in the rat on inotropic stimulation of the adult heart, Cardiovasc Res. 65, 861-8.

12. Pei, J. M., Yu, X. C., Fung, M. L., Zhou, J. J., Cheung, C. S., Wong, N. S., Leung, M. P. & Wong, T. M. (2000) Impaired G(s)alpha and adenylyl cyclase cause beta-adrenoceptor desensitization in

chronically hypoxic rat hearts, Am J Physiol Cell Physiol. 279, C1455-63.

13. Zhu, W., Zeng, X., Zheng, M. & Xiao, R. P. (2005) The enigma of beta2-adrenergic receptor Gi signaling in the heart: the good, the bad, and the ugly, Circ Res. 97, 507-9.

14. Liu, R., Ramani, B., Soto, D., De Arcangelis, V. & Xiang, Y. (2009) Agonist dose-dependent phosphorylation by protein kinase A and G protein-coupled receptor kinase regulates beta2 adrenoceptor coupling to G(i) proteins in cardiomyocytes, J Biol Chem. 284, 32279-87.

15. Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmuller, U. & Timmer, J. (2009) Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood, Bioinformatics. 25, 1923-9.

16. Ahlquist, R. P. (1948) A study of the adrenotropic receptors, Am J Physiol. 153, 586-600. 17. Toews, M. L., Harden, T. K. & Perkins, J. P. (1983) High-affinity binding of agonists to beta-adrenergic receptors on intact cells, Proc Natl Acad Sci U S A. 80, 3553-7.

18. Green, S. A., Holt, B. D. & Liggett, S. B. (1992) Beta 1- and beta 2-adrenergic receptors display subtype-selective coupling to Gs, Mol Pharmacol. 41, 889-93.

19. Cedersund, G. (2012) Conclusions via unique predictions obtained despite unidentifiability--new definitions and a general method, FEBS J. 279, 3513-27.

20. Brännmark, C., Palmér, R., Glad, S. T., Cedersund, G. & Strålfors, P. (2010) Mass and information feedbacks through receptor endocytosis govern insulin signaling as revealed using a parameter-free modeling framework, J Biol Chem. 285, 20171-9.

(20)

21. Nyman, E., Brännmark, C., Palmér, R., Brugård, J., Nyström, F. H., Strålfors, P. & Cedersund, G. (2011) A hierarchical whole-body modeling approach elucidates the link between in Vitro insulin signaling and in Vivo glucose homeostasis, J Biol Chem. 286, 26028-41.

22. Tafintseva, V., Tondel, K., Ponosov, A. & Martens, H. (2014) Global structure of sloppiness in a nonlinear model, Journal of Chemometrics. 28, 645-655.

23. Eisenberg, M. C. & Hayashi, M. A. (2014) Determining identifiable parameter combinations using subset profiling, Math Biosci. 256, 116-26.

24. Saucerman, J. J., Brunton, L. L., Michailova, A. P. & McCulloch, A. D. (2003) Modeling beta-adrenergic control of cardiac myocyte contractility in silico, J Biol Chem. 278, 47997-8003.

25. Saucerman, J. J. & McCulloch, A. D. (2004) Mechanistic systems models of cell signaling networks: a case study of myocyte adrenergic regulation, Prog Biophys Mol Biol. 85, 261-78.

26. Saucerman, J. J., Zhang, J., Martin, J. C., Peng, L. X., Stenbit, A. E., Tsien, R. Y. & McCulloch, A. D. (2006) Systems analysis of PKA-mediated phosphorylation gradients in live cardiac myocytes, Proc

Natl Acad Sci U S A. 103, 12923-8.

27. Yang, J. H. & Saucerman, J. J. (2012) Phospholemman is a negative feed-forward regulator of Ca2+ in beta-adrenergic signaling, accelerating beta-adrenergic inotropy, J Mol Cell Cardiol. 52, 1048-55. 28. Holland, D. O., Krainak, N. C. & Saucerman, J. J. (2011) Graphical approach to model reduction for nonlinear biochemical networks, PLoS One. 6, e23795.

29. Bondarenko, V. E. (2014) A compartmentalized mathematical model of the beta1-adrenergic signaling system in mouse ventricular myocytes, PLoS One. 9, e89113.

30. Xin, W., Tran, T. M., Richter, W., Clark, R. B. & Rich, T. C. (2008) Roles of GRK and PDE4 activities in the regulation of beta2 adrenergic signaling, J Gen Physiol. 131, 349-64.

31. Shankaran, H., Wiley, H. S. & Resat, H. (2007) Receptor downregulation and desensitization enhance the information processing ability of signalling receptors, BMC Syst Biol. 1, 48.

32. Nyman, E., Fagerholm, S., Jullesson, D., Stralfors, P. & Cedersund, G. (2012) Mechanistic

explanations for counter-intuitive phosphorylation dynamics of the insulin receptor and insulin receptor substrate-1 in response to insulin in murine adipocytes, FEBS J. 279, 987-99.

33. Brännmark, C., Nyman, E., Fagerholm, S., Bergenholm, L., Ekstrand, E. M., Cedersund, G. & Strålfors, P. (2013) Insulin signaling in type 2 diabetes: experimental and modeling analyses reveal mechanisms of insulin resistance in human adipocytes, J Biol Chem. 288, 9867-80.

34. Keely, S. L. & Corbin, J. D. (1977) Involvement of cAMP-dependent protein kinase in the regulation of heart contractile force, Am J Physiol. 233, H269-75.

35. Devic, E., Xiang, Y., Gould, D. & Kobilka, B. (2001) Beta-adrenergic receptor subtype-specific signaling in cardiac myocytes from beta(1) and beta(2) adrenoceptor knockout mice, Mol Pharmacol. 60, 577-83.

36. Fisher, G. W., Adler, S. A., Fuhrman, M. H., Waggoner, A. S., Bruchez, M. P. & Jarvik, J. W. (2010) Detection and quantification of beta2AR internalization in living cells using FAP-based biosensor

technology, J Biomol Screen. 15, 703-9.

37. Hodding, G. C., Jann, M. & Ackerman, I. P. (1980) Drug withdrawal syndromes-- a literature review,

West J Med. 133, 383-91.

38. Schmidt, H. & Jirstrand, M. (2006) Systems Biology Toolbox for MATLAB: a computational platform for research in systems biology, Bioinformatics. 22, 514-5.

39. Schmidt, H. (2007) SBaddon: high performance simulation for the Systems Biology Toolbox for MATLAB, Bioinformatics. 23, 646-7.

40. Klipp, E., Liebermeister, W., Wierling, C., Kowald, A., Lehrach, H. & Herwig, R. (2009) Systems

Biology, Wiley-Blackwell.

41. Cedersund, G. & Roll, J. (2009) Systems biology: model based evaluation and comparison of potential explanations for given biological data, FEBS J. 276, 903-22.

42. Kreutz, C., Raue, A. & Timmer, J. (2012) Likelihood based observability analysis and confidence intervals for predictions of dynamic models, BMC Syst Biol. 6, 120.

(21)

(22)

Figure legends

Figure 1 - Chicken ventricular tissue displays decreased force of contraction at high concentrations of adrenaline (ADR).

A The experimental protocol for the ADR doses (0.01 – 10000 nM) and corresponding times of addition.

B The resulting dose-response curve of force of contraction in cardiac tissue normalized to the basal response (mean±SEM, n=6). Indicated with brown color are the last measurements before a new dose is added. These measurements are used in Figure 1D as indicated with dashed lines.

C Individual dose-response experiments with corresponding data in Figure 1B. Indicated in yellow and orange are 2 curves that do not have an obvious decrease in response at high doses of ADR.

D Sigmoidal curves with variable slope were fitted to the data (Materials and Methods). Dashed lines indicate the transition of data from Figure 1B.

E The relation between the EC50 value of the sigmoidal curves in Figure 1D and their respective agreement with data gives a 95 % confidence interval for EC50 of 46-191 nM. The statistical threshold χ2(7-4,0.05)=7.8 is indicated with a red line.

Figure 2 - A minimal mathematical model of beta-adrenergic signaling.

A A sketch of the minimal mathematical model that is used in the analysis of data. H represents the agonist, β represents the basal state of the inactive receptor, βact represents

(23)

represents the amount of cyclic AMP which herein is directly translated to the force of contraction. The rate parameters k1-k4, kbasal, kH, and cAMP0 are indicated near the corresponding rate arrows.

B The agreement between model simulations with an approximation of all found acceptable parameters (Materials and Methods) and the experimentally determined dose-response curve (from Figure 1B). The best agreement is indicated in bold. Data

(mean±SEM) are shown with asterisks and model simulations with lines.

Figure 3 - A dynamic overshoot behavior is predicted as the underlying cause to the decrease in the dose-response curve.

A Time-resolved simulations of the minimal model with a step input of different

adrenaline (ADR) doses at time = 5 minutes. The different lines represent the different input doses of 0.01 – 10000 nM ADR.

B Accumulative time-resolved model simulations of the minimal model with continuous responses to ADR. The ADR concentration in the simulation is increased according to the experimental protocol in Figure 1A, but for each dose the simulation continues to time = 22 minutes. The different lines represent the different concentrations (0.01 – 10000 nM ADR). Indicated in bold is the resulting dose-response curve with a decreasing response for high doses.

C The minimal model predicts an overshoot behavior for a single dose of 1000 nM added ADR at time = 5 minutes. The different lines represent simulations for an approximation of all found acceptable parameters (Materials and Methods) and accounts for the uncertainty of the prediction.

(24)

D Model simulations of the steady-state time-resolved dose-response curve for ADR stimulated force of contraction. The different lines represent simulations for an

approximation of all found acceptable parameters and accounts for the uncertainty of the prediction.

E The corresponding sigmoidal dose-response curve gives EC50 of 44-369 nM ADR for all acceptable parameter sets. The different lines represent simulations for an approximation of all found acceptable parameters and accounts for the uncertainty of the prediction.

Figure 4 - Time-resolved isoproterenol (ISO) stimulated cardiac tissue confirm the overshoot behavior.

A Time-resolved response data for a serial addition of 100 nM ISO for 10 minutes after 5 minutes of rest, washing for 5 minutes, and 100 nM ISO for another 10 minutes

(mean±SEM, n=11).

B The fit between the model and the time-resolved data is poor (red line), mainly due to the first data point after first addition of ISO (marked with a red circle). The introduction of a time-delay parameter gives a statistically good enough fit with data (green line)

(χ2

(minimal model)=18.7 < 31). Data (mean±SEM) are shown with blue asterisks and model simulations with red and green lines.

C Profile likelihood analysis of the time-delay parameter gives 0.65 minutes = 39 seconds as the optimal value of the parameter.

Figure 5 - Accounting for the difference between the two agonists isoproterenol (ISO) and adrenaline (ADR).

(25)

Profile-likelihood analysis of the relation between agreement with data and the ratio between the association constant of ISO and ADR (kratio). The best agreement with both ISO and ADR data is achieved for kratio = 5.5. The interval of acceptable values is 3<kratio<9.

Figure 6 - Recalculating a model-based EC50 confidence interval for adrenaline (ADR)

stimulated chicken ventricular tissue.

A Model simulations in agreement with the dose-response and the time-resolved data. The different lines represent simulations of the model with different acceptable parameter sets. Data (mean±SEM) are shown with asterisks and model simulations with lines.

B Model simulations of a steady-state experiment with 30 minutes waiting time before new additions of ADR. The different lines are simulations with the different acceptable parameter sets in Figure 6A.

C Model simulations of the steady-state dose-response curve for the acceptable parameter sets in Figure 6A. The resulting calculated EC50 = 193-493 nM ADR for all found

acceptable parameter sets.

D A 95 % confidence interval for the model-based EC50 gives values of 180-525 nM under the χ2 threshold of 95 (indicated with a black line). Specific searches in the space of

parameters were performed to obtain the confidence interval (Materials and Methods).

E Slightly overlapping but substantially different confidence interval from the classical interpretation of EC50 in Figure 1E (red) and the model-based interpretation of EC50 in Figure 6D (black).

(26)

0 5 10 15 20 25 1 1.2 1.4 1.6 1.8 Time (minutes) Force of contraction

, fold over basal

Experimental data 10−2 10−1 100 101 102 103 104 1 1.2 1.4 1.6

Force of contraction, fold over basal

0 5 10 15 20 25 1 1.5 2 Time (minutes) Force of contraction

, fold over basal

Individual experiments 0 5 10 15 20 25 0 10−2 100 102 104 Time (minutes) Doses of adrenaline Adrenaline (nM)

Figure 1

A

C

B

D

E

Agreemant with data (cost) 4 5 6 7 8 9 102 103 101 Classical method 46-191 nM Classical dose-response

(27)

Figure 2

A

B

β

H

βact

βphos

cAMP

kH kbasal k4 k2 k1 k3 cAMP0 0 10 20 1 1.5 Time, minutes

Simulation of acceptable parameters

(28)

0 20 40 60 80 1

1.5

Time, minutes

Single dose simulation (1000 nM)

0 200 400 600

1 1.25

Time, minutes

Steady-state simulations 0 10 20 1 1.5 Time, minutes

Single dose stimulations (0.01-1000 nM) 0 10 20 1 1.5 Time, minutes Force of contraction,

fold over basal

Accumulative stimulations (0.01-10000 nM) −2 −1 0 1 2 3 4 0 50 100 Response, % of max Dose-response curve

Figure 3

A

B

C

D

E

(29)

0 10 20 30 0.8 1 1.2 1.4 1.6 1.8 Time (minutes)

Time delay=0, cost=44.6 Time delay=0.65, cost=18.7

0.2 0.4 0.6 0.8 1 18 19 20 21 22 23

Time delay (minutes)

Agreement with data (cost)

0 10 20 30 0.8 1 1.2 1.4 1.6 Time (minutes)

ISO, 100 nM ISO, 100 nM Experimental data

Figure 4

A

B

(30)

2 4 6 8 10 88 89 90 91 92 93 94

parameter kratio (dimensionless)

Agreement with data (cost)

(31)

0 10 20 1

1.5

Time (minutes)

Force of contraction fold over basal

0 10 20

1 1.5

Time (minutes)

0 200 400 600

1 1.5

Time (minutes)

Steady-state simulations 100 −2 10−1 100 101 102 103 104 50 100 Adrenaline (nM) Response, % of max Dose-response curve 101 102 103 88 90 92 94 96 EC50 (nM)

Agreement with data (cost) 180-525 nM Model-based method 101 102 103 90 94 EC50 (nM) cost cost 180-525 Classical method

Figure 6

A

B

C

D

E

5 6 7 8 Model-based method 46-191

(32)

Table I. The found acceptable values of the model parameters compared to the chosen allowed range in the analysis of the final model.

Name of parameter Best value found Minimal acceptable value found Maximal acceptable value found Minimal value allowed Maximal value allowed kbasal 0.112 0.0719 0.228 0.0001 10000 kH 0.000775 0.000556 0.00108 0.0001 10000 k1 0.116 0.0842 0.153 0.05 10000 k2 0.837 0.515 1.38 0.0001 10000 k3 1.23 0.0001 1303 0.0001 10000 k4 218 5.59 10000 0.0001 10000 cAMP0 18.3 0.000667 10000 0.0001 10000 kratio 5.55 3.86 8.31 3 9 kscale 1 1 1 1 1 time delay 0.65 0.65 0.65 0.65 0.65