The Application of Dynamic Models to the Exploration of β1-AR Overactivation as a Cause of Heart Failure

High titer of β1-adrenoreceptor autoantibodies (β1-AA) has been reported to appear in heart failure patients. It induces sustained β1-adrenergic receptor (β1-AR) activation which leads to heart failure (HF), but the mechanism is as yet unclear. In order to investigate the mechanisms causing β1-AR non-desensitization, we studied the beating frequency of the neonatal rat cardiomyocytes (NRCMs) under different conditions (an injection of isoprenaline (ISO) for one group and β1-AA for the other) and established three dynamic models in order to best describe the true relationships shown in medical experiments; one model used a control group of healthy rats; then in HF rats one focused on conformation changes in β1-AR; the other examined interaction between β1-AR and β2-adrenergic receptors (β2-AR). Comparing the experimental data and corresponding Akaike information criterion (AIC) values, we concluded that the interaction model was the most likely mechanism. We used mathematical methods to explore the mechanism for the development of heart failure and to find potential targets for prevention and treatment. The aim of the paper was to provide a strong theoretical basis for the clinical development of personalized treatment programs. We also carried out sensitivity analysis of the initial concentration β1-AA and found that they had a noticeable effect on the fitting results.


Introduction
The incidence of heart failure (HF) is increasing year by year throughout the world, as is the cost of treating it [1]. Accurately recognising the cardiac signals associated with HF not only would prevent the progress of chronic HF, but also can guide individual treatment and positively influence the normal progression of the disease [2]. A variety of mechanisms are involved in HF progression, including nervous system inhibition. Heart failure itself causes cardiovascular irregularity due to material imbalance and cell damage [3]. Studies have shown that sustained activation of the sympathetic nervous system (SNS) is the core mechanism trigger of heart failure [4]. The main way in which the overactivation of the SNS occurs is the overactivation of the cardiomyocytes at the 1 -adrenergic receptor ( 1 -AR) [5]. 1 -adrenergic receptor autoantibodies ( 1 -AA) were found in the serum of patients with dilated cardiomyopathy in a study [6]. As 1 -AA can bind with and activate 1 -AR, it has a 1 -AR agonistlike effect [7]. Excessive activation of 1 -AR leads to impaired cardiac function [8]. A study found that, in addition to catecholamine 1 -AR agonist ISO, 1 -adrenoceptor ( 1 -AA) autoantibodies that could cause continuous activation of 1 -AR [9] were detected in 40%-60% of HF patients [8]. 1 -AR and 2 -AR ( 1 -adrenergic receptors) link through the Cterminal, which are within cardiac cells, and belong to the G protein-coupled receptor family [10]. 1 -AR and 2 -AR can form heterodimers [11], where 2 -AR downstream couples with stimulatory G protein ( ) and inhibitory G protein ( ) [12]. Another study [13] found that, in the early stages 2 Computational and Mathematical Methods in Medicine of heart failure, when 1 -AR is activated, it can cause the 2 -AR downstream signal to be converted from to , and is activated by the activity of G protein-coupled receptor kinase 2 (GRK2), which then causes 1 -AR phosphorylation, as well as -arrestin binding to form endocytosis, after which 1 -AR endocytosis [13] inhibits the overactivation of 1 -AA. There are two main mechanisms involved in the continuous activation of 1 -AR: 1 -AR conformation changes and interaction between 1 -AR and 2 -AR. However, the mechanisms for how 1 -AA induces continuous activation of 1 -AR are not entirely clear.
To identify the core molecular mechanism of 1 -AR nondesensitization, we adopted a novel approach by establishing differential dynamic models. We used ++ to estimate the parameters, obtained modified 60-minute snapshots highlighting the experimental data, and made predictions on how the images would progress over the next 60 minutes. Then we took advantage of AIC to select the optimal model of 1 -AR non-desensitization induced by 1 -AA, and the corresponding mechanism. 2 -AR plays an important role in the activation of 1 -AR, which can be used to estimate cardiac function, prevent deterioration of cardiac function, and improve the quality of life of the patients with HF.
This paper is organized as follows. In Section 2, we established the models and listed the differential equations separately. Section 3 focused on the analysis of the models and determination of the parameters. We select the optimal model to identify the key molecular mechanism, by applying the Akaike information criterion (AIC). Results on the sensitivity analysis of parameters are given in Section 4. Section 5 is the conclusion. Figure 1 shows the interrelationships between the various substances. It reflects the receptor desensitization when the 1 -AR is combined with ISO and non-desensitization when the 1 -AR is combined with 1 -AA, causing conformational changes and interaction between the 1 -AR and the 2 -AR, eventually leading to HF.

Extraction and Detection of Neonatal Rat Cardiomyocytes.
Laboratory animal medicine of Capital Medical University provided 20 male newborn rats born 0-3 days as raw materials for extracting neonatal rat cardiomyocytes. According to previous method [14], the details of myocardial cell isolation and culture of neonatal rats are as follows: (1) open the sternum, expose the heart, clip it with tweezers, and wash in cold phosphate-buffered saline (PBS); (2) remove excess connective tissue from the washed heart and cut with ophthalmic scissors; (3) the heart fragments were aspirated into a centrifuge tube, and cold PBS was added and then centrifuged at 1000 for 5 ; (4) remove the centrifuge tube and vacuum suction PBS; (5) add 1 of 0.25 % trypsin and 1 of 0.25 % collagenase, blow vigorously, set in a 37 ∘ C water bath, and shake for 20 ; (6) add 1 fetal bovine serum to stop digestion; (7) digestion was collected by adding 2 of DMEM containing 10% fetal bovine serum (FBS), 1000 , and after centrifugation for 10 , the supernatant was discarded and fresh medium was added for differential adherence; (8) after adherence for 1 ℎ, the culture medium was transferred to a new centrifuge tube and centrifuged at 1000 for 5 . Cells were collected and the supernatant was discarded. Fresh medium was added and transferred to a 6-well plate for cell culture and 2 of DMEM lowglucose medium was added to each well. Detection of the beating frequency of NRCMs [15] was as follows: the culture medium was replaced on the day of the experiment and stably incubating in a 37 ∘ C cell culture incubator for 30 minutes, the 6-well plate was placed on a constant-temperature table of inverted microscopy, and a total of 10 fields of view of 3 sixwell plates were randomly observed. Each field was measured for 30 seconds at a time, and the number of synchronized contractions of an isolated single cell or a group of cardiomyocytes in the untreated group was measured.
Then, we used 0.1 of ISO, 1 -AA, and IgG (immunoglobulin G) (eliminating the effects of 1 -AA itself) to stimulate, respectively. After that, the beating frequency of NRCMs was measured by Live Cell Imaging System offered by Medical Sciences Center Lab of Capital Medical University. Specifically, we measured the frequency in beat per minute of the cardiomyocytes after stimulation by live cell workstation and visualization at 63-fold magnification. Before being exposed to drugs, the cells had been stabilized for 10 in the system. The present study complies with the recommendations in the Guide for the Care and Use of Laboratory Animals protocol, NIH guidelines (Guide for the Care and Use of Laboratory Animals), and conformed to AVMA Guidelines on Euthanasia.
Tables 1 and 2 showed the measured data of beating frequency of NRCMs added to ISO and 1 -AA. We used the beat frequency of NRCMs at time 0 as the base value, which was recorded as 0. When the standard deviation was greater than or equal to 4, the data was appropriately adjusted within the scope of standard deviation. Figure 2 reflects the positive correlation between concentration and beating frequency of NRCMs. Also, Martinsson et al. [16] studied the relationship between ISO and heart rate in their experiment results.

The Models.
There are many molecular mechanisms causing 1 -AR non-desensitization, among which conformation changes and interaction are the most likely ones. In order to clarify the most likely molecular mechanism, we established dynamical models including control model and experiment models to study the specific molecular mechanisms.

Explanation of Interaction
Definition 1 (a protein-to-protein interaction ( ) [17]). Proteins rarely act alone as their functions tend to be regulated. Many molecular processes within a cell are carried out by molecular machines that are built from a large number of protein components organized by their (protein-protein interactions). In the experimental groups, 1 -AR and 2 -AR belong to the G protein-coupled receptor family that are connected through the C-terminus [10] which belongs to PPI. Although 1 -AA does not directly bind to 2 -AR, 1 -AA will indirectly affect the conformation of 2 -AR by activating 1 -AR based on the studies of laboratory animal medicine of Capital Medical University. Therefore, 1 -AA can interfere with their dimerization which reflects the fact that 2 -AR can affect the persistence of 1 -AR activation. Figure 3: The block diagram of control model. ISO (L) activates 1 -AR ( 1 ) receptor generating intermediate complexes 1 causing endocytosis of 1 -AR, so that it can no longer come into contact with the protected ligand; then the complexes LB produce a new substance P and structurally changed receptor 1 -AR ( 1 ) that returns to the cell surface to terminate the signal. Figure 4: The block diagram of conformation changes model. 1 -AA (I) activates 1 -AR ( 1 ) receptor generating intermediate complexes 1 with conformation changes. The production of 1 continues to activate the signal, and then it produces a new substance P and structurally changed receptor 1 -AR ( 1 ).

The Establishment of the Model Block Diagrams.
The possible molecular mechanisms about 1 -AR overactivation are conformation changes and interaction. To determine a more specific molecular mechanism, we propose the dynamical models including control model and experiment models, where experiment models are composed of conformation changes model and interaction model. In order to facilitate the study and simplify the reaction diagram, we use corresponding letters instead of reactants and the definition of the corresponding parameters in dynamical models; see Table 3.
(1) The Block Diagram of Control Model (See Figure 3). In Figure 3, 1 represents the combined velocity of L and 1 , −1 represents the reverse reaction velocity, 2 represents the velocity of 1 decomposition, −2 is the reverse reaction velocity, and 1 and 2 represent the degradation velocity of 1 and P, respectively. 3 represents the velocity at which undegraded 1 returns to the cell surface and participates in the reaction again. Figure 4). In Figure 4, 1 represents the combined velocity of I and 1 , −1 represents the reverse reaction velocity, 2 represents the velocity of 1 decomposition, −2 is the reverse reaction velocity, and 1 and 2 represent the degradation velocity of 1 and P, respectively. 3 represents the velocity at which undegraded 1 returns to the cell surface and participates in the reaction again. Figure 5). In Figure 5, 1 represents the combined velocity of I and 1 2 , −1 represents the reverse reaction velocity, 2 represents the velocity of 1 2 decomposition, 3 represents the velocity of ( 1 ) decomposition, and 1 and 2 represent the degradation velocity of 1 and P, respectively. 3 represents the velocity at which undegraded 1 returns to the cell surface and participates in the reaction again. Also, 1 -AR and 2 -AR  Figure 5: The block diagram of interaction model. We also take 2 -AR ( 2 ) into consideration in the reaction based on medical experiments. According to the studies of laboratory animal medicine of Capital Medical University, they found that 1 -AA had a unique characteristic that can inhibit heterodimerization of 1 / 2 -AR and thus lead to decompose of 1 -AR and 2 -AR that are originally connected, resulting in endocytosis inhibition and sustained activation of the signal. Therefore, there is no conjugation of 1 -AR and 2 -AR ( 1 2 ) in the products. Besides, 1 -AA (I) cannot directly bind to 2 -AR ( 2 ), so there is no 2 in the products. Thus, we consider that 1 2 (I and 2 are combined to 1 ) are intermediate complexes and 1 , 2 are the products with structure change.

(3) The Block Diagram of Interaction Model (See
belong to the G protein-coupled receptor family [10], and thus we assume that they have the same degradation velocity 1 .
Next, we established three mathematical models as shown in Figures 3, 4, and 5 and established the corresponding ordinary differential equations based on the relevant theoretical knowledge of biochemical reaction models in cell and molecular biology [18,19].

Corresponding Differential Equations
(1) The Control Model. we propose a control model for healthy rats based on Figure 3. The time evolution of control model is described by five coupled differential equations.
Computational and Mathematical Methods in Medicine 5  (2) The Conformation Changes Model. We have known that 1 -AA and ISO compete for different binding sites of 1 -AR, which belongs to the competitive inhibition in different positions of 1 -AR [10]. Assuming that there is no reaction between 1 -AA and 2 -AR, the time evolution of conformation changes model based on Figure 4 is described by five coupled differential equations.
where the initial conditions of the conformation changes , Next, we will use second order Runge-Kutta method to estimate the parameters of systems (1), (2), and (3) and then select the optimal model of 1 -AR non-desensitization.

Parameter Fitting and Graphs Analysis. Second order
Runge-Kutta method is used to solve the ODE systems obtained in our models. Our parameter fitting is an optimization process, which minimizes the squared summation of . We also denote the total squared summation by . Obviously is a function of all the parameters that require fitting. Given a set of parameters, second order Runge-Kutta method gives the value of [ ] and thus . We use steepest descent method to find the minimum of . In this process, the gradient of over all parameters is calculated numerically. For example, for positive parameter , In our program Δ = 1.0 − 6, and the time step in the Runge-Kutta method is chosen to be 0.01 min.

Computational and Mathematical Methods in Medicine
We use the measured values in Table 1 to do the parameter fitting for the control model and the values in Table 2 for the experiment models. Note that the measured data in the tables is not the exact concentration of [ ]. Instead, we assumed that the measured data is proportional to the concentration of [ ]. In other words, we introduce an extra parameter , and [ ] = , where is the measured data in the tables. In the optimization process, we trait in the same way we trait other parameters. And is also get fitted. Combining Figure 2 and parameter fitting method, we determined that the concentration is directly proportional to the heart rate of NRCMs and the proportionality factor is 0.0025.
Next, we analyze the fitting of the graphs and the prediction of the different time points, as shown in Figures 6-8. Figure 6(a) shows the fitting results of control group model from 0 to 60 , in which the points are experimentally measured and the curve is the solution of dynamical control group model of [ 1 ]. And the fitting curve finally declined which reflects the fact that 1 -AR desensitizes when reacting with ISO in healthy rats. However, the situation was quite different with HF rats. Figure 6(b) displays the prediction result from 60 to 120 . It predicts that [ 1 ] continues to decrease, eventually tends to zero, and becomes stable, which clearly showed the mechanism of 1 -AR desensitization. Figure 7(a) presents the fitting results of the conformation changes model from 0 to 60 , in which the points are experimentally measured and the curve is the solution of dynamical conformation changes model of [ 1 ]. Figure 7(b) is the prediction results from 60 to 120 . The [ 1 ] curve eventually stabilized in (b), which clearly showed the mechanism of 1 -AR non-desensitization. Figure 8(a) presents the fitting results of the interaction model from 0 to 60 , in which the points are experimentally measured, the curve is the solution of dynamical interaction model of [ 1 2 ], and the fitting curve finally stabilized. Figure 8(b) is the prediction results from 60 to 120 . The upward trend of the [ 1 2 ] curve reflects the effects of 1 -AA and 2 -AR on the non-desensitization of the 1 -AR.
Then, we analyze the reaction speeds and Table 4 lists the reaction velocities of three models. Since the units of positive reaction velocities and reverse reaction velocities are different, they cannot be directly compared. Therefore, we simulated the reaction speeds figures of the control group and the experimental groups, and the reaction speed is the product of the reaction velocity and concentration. See Figure 9.
As can be seen from Figure 9, the positive and negative reaction speeds of the control model and experimental models are greater than those of the reverse reaction, and the positive reaction speeds decrease rapidly. In the course of the decrease of the positive reaction speeds, the speed of reversed reaction continues to rise. Finally, their tendency gradually slows down. [20][21][22] (AIC) is a measure of the goodness of fit of an estimated statistical model and a tool for model selection. Given a data set, several competing models can be ranked according to their AIC, with the one with the lowest AIC being the best. The AIC is defined as follows:

Model Selection. Akaike's information criterion
where RSS is the residual sum of squares, is the amount of the real data, and is the total amount of estimated parameters in the model. AIC is AIC with a second order correction for small sample sizes ( / < 40), to start with Moreover, AIC difference is defined as Δ = − , where is the AIC or AIC of the th model, and is the AIC or AIC of the model with minimal AIC. If 0 ≤ Δ ≤ 2, the th model has substantial data support. If 4 ≤ Δ ≤ 7, the th model has considerable less data support. If Δ ≥ 10, the th model has essentially no data support.
Since the / = 10/4 < 40, we use AIC to select the model, which determines the mechanism of the 1 -AR nondesensitization. In order to verify the validity of the models, we apply AIC to view the residuals between theoretic and real data under three models. We calculated the AIC values of the conformation model and interaction model which are, respectively, 1 ≈ −7.287, 2 ≈ −8.494. By comparing the AIC values calculated in the two models above, we conclude easily that 1 > 2 , ≈ −8.494. It is obvious that the possibility of interaction is the largest, which conforms to the analysis of graphs and parameter.

Sensitivity Analysis
We determined that the interaction model between 1 -AR and 2 -AR is the optimal model for heart failure through the AIC criterion. Next, when the initial concentration of 1 -AA is different, we analyzed the figures and data (measured by Capital Medical University) fitting of the interaction model between 1 -AR and 2 -AR. From Figure 10, we can see that the initial concentration of 1 -AA has a certain effect on the fitting result. Table 5 shows the reaction velocity for the interaction model at different concentrations. Figure 10, respectively, showed the graphs of [ 1 2 ] when the initial concentrations of 1 -AA are 1 × 10 −9 / , 1 × 10 −8 / , and 1 × 10 −6 / in 10 minutes, while the concentration of 1 2 remains unchanged at 1.8×10 −9 / . And we can see that when the concentration of 1 -AA is 1 × 10 −9 / , the fitting result is the worst, followed by 1 × 10 −6 / and 1 × 10 −8 / . As can be seen from Figure 8, when concentration of 1 -AA is 1 × 10 −7 / , the fitting result is the best. Therefore, the optimum concentration of mice monoclonal 1 -AA leading to the NRCMs beating frequency was 0.1 (detected from 1 to 1 ).
Computational and Mathematical Methods in Medicine

Conclusion
We analyzed the two types of mechanisms that cause HF; it can be concluded that conformation changes and interactions between 1 -AR and 2 -AR both contribute to the HF progression in rats and the interaction has the greater impact on it. In our study, for the first time, the model block diagrams described the causation of HF. By comparing the experimental group models with the control group model and then calculating the AIC , we concluded that the interaction model between 1 -AR and 2 -AR was the optimal model, which leads to the important role of 2 -AR in the progression of HF. 2 -AR and 1 -AR connect together and participate in a reaction, which then causes conformation changes in receptor 1 -AR. This results in endocytosis inhibition, which ultimately leads to 1 -AR non-desensitization and HF. Thus, 2 -AR is an important consideration in guiding the treatment of HF.
Previous studies have found that 1 -AR and 2 -AR could form heterodimers, although they also found that the autoantibody 1 -AA did not bind to 2 -AR [15]. We found that 1 -AA could inhibit heterodimerization of 1 / 2 -AR, observing laboratory experiments on animals at Capital Medical University. Therefore, we speculated that 1 -AA might indirectly affect the 2 -AR-C-terminal conformation and interfere with heterodimerization of 1 / 2 -AR. This would lead to 1 -AR endocytosis inhibition, to continuous signal activation, and ultimately to heart failure.
Because the interaction between 1 -AR and 2 -AR was found to be the trigger mechanism of heart failure progression, receptor 2 -AR should be further studied to determine if it could play a role in heart failure treatment. In the presence of 1 -AA-positive heart failure, 2 -AR can be used as a drug and therapeutic target to improve the interaction between 1 -AR and 2 -AR, thereby decreasing the risk of further heart failure.

Data Availability
The study conformed to AVMA Guidelines on Euthanasia and the Guide for the Care and Use of Laboratory Animals protocol published by the Ministry of Education of People's Republic of China. All studies were approved by Capital Medical University Committee on Animal Care. According to the data provided by Capital Medical University, we extracted the original data and established the relationship between experimental data and mathematical models to facilitate our research.

Ethical Approval
Our research was approved by the Institutional Animal Care and Use Committee of the Capital Medical University. The investigators understand the ethical principles under which the journal operates, and the study complies with the journal's animal ethics checklist (ethical number is AEEI-2016-053).

Conflicts of Interest
The authors declare that there are no conflicts of interest.
Talent Introduction Research Fund of Taiyuan University of Technology (tyut-rc201317a), and China Scholarship Council.

Supplementary Materials
Supplementary 1. The picture in the supplementary material reflects the increased beating frequency of NRCMs stimulated with different concentrations of 1 -AA. As can be seen, 0.1 1 -AA is the optimal concentration for the increase of beating frequency in NRCMs, which is consistent with the results of the sensitivity analysis.
Supplementary 2. The video file in the supplementary material mainly indicates that we successfully extracted NRCMs from the experiments and measured the beating frequency.