Multiscale Modelling of -Adrenergic Stimulation in Cardiac Electromechanical Function

β-adrenergic receptor stimulation (β-ARS) is a physiological mechanism that regulates cardiovascular function under stress conditions or physical exercise. Triggered during the so-called “fight-or-flight” response, the activation of the β-adrenergic receptors located on the cardiomyocyte membrane initiates a phosphorylation cascade of multiple ion channel targets that regulate both cellular excitability and recovery and of different proteins involved in intracellular calcium handling. As a result, β-ARS impacts both the electrophysiological and the mechanical response of the cardiomyocyte. β-ARS also plays a crucial role in several cardiac pathologies, greatly modifying cardiac output and potentially causing arrhythmogenic events. Mathematical patient-specific models are nowadays envisioned as an important tool for the personalised study of cardiac disease, the design of tailored treatments, or to inform risk assessment. Despite that, only a reduced number of computational studies of heart disease have incorporated β-ARS modelling. In this review, we describe the main existing multiscale frameworks to equip cellular models of cardiac electrophysiology with a β-ARS response. We also outline various applications of these multiscale frameworks in the study of cardiac pathology. We end with a discussion of the main current limitations and the future steps that need to be taken to adapt these models to a clinical environment and to incorporate them in organ-level simulations.


Introduction
β-adrenergic receptor stimulation (β-ARS) is a physiological response mechanism that plays a fundamental role in the regulation of cardiomyocyte activity, producing a positive inotropic (enhanced contraction), lusitropic (faster relaxation), and chronotropic (increased heart rate) effect. Such a multifactorial response is triggered via the activation of the β-adrenergic receptors by the sympathetic nervous system, under either stress conditions or physical exercise, and is also known as the "fight-or-flight" response.
β-adrenergic receptors were first described by Lands et al. in the late 1960s [1,2]. They are situated on the cardiomyocyte membrane and react to different neurotransmitters (norepinephrine, epinephrine) and drugs (isoprenaline). When binding with the adrenergic receptor takes place, it starts a reaction cascade where different cellular substrates become phosphorylated, affecting their individual roles in the overall excitation-contraction coupling. As a result, under healthy physiological conditions, the cardiac action potential shortens, while the intracellular calcium transient exhibits an increased amplitude and a faster decay rate as the main manifestations of β-ARS at the cellular level [3,4]. However, the large number of components involved in the β-adrenergic cascade, and the complexity of these subcellular processes and interactions, make β-ARS signalling highly sensitive to cellular changes and to pathological perturbations. As a result, β-ARS plays a main role in a considerable number of heart diseases [5], and it is well established as an important contributor to cardiomyocyte arrhythmogenicity [6][7][8].
In particular, the β-ARS response has direct effects on the ion channels and pumps of the cell membrane (and, therefore, on intracellular ionic concentrations) regulating calcium intake, intracellular calcium handling, calcium extrusion, and cellular repolarisation. Impairments in the balance between these carefully orchestrated processes can affect heart function and render its constituent cardiomyocytes susceptible to proarrhythmic events, such as early and delayed afterdepolarisations (EADs and DADs, respectively). Such afterdepolarisations under β-ARS are common proarrhythmic manifestations in isolated cardiomyocytes from patients of different pathologies, especially those characterised by action potential and calcium transient abnormalities (such as hypertrophic cardiomyopathy [9], long QT syndrome [10,11], or catecholaminergic polymorphic ventricular tachycardia [12]). An overstimulated β-ARS response is also one of the main mechanisms of cardiac hypertrophy, coronary artery disease, or stroke events [13]. The overexpression of the β-adrenergic response has also been linked to the onset of cardiac hypertrophy or the generation of fibrotic tissue [14,15]. The appearance of these structural changes can lead to the creation of re-entry pathways in myocardial tissue, which also contribute to the generation of self-sustained arrhythmias. Induced arrhythmias are also a common manifestation in heart failure. In chronic heart failure, cardiac remodelling at the structural level can affect the pathways involved in the β-ARS response [16]. As a result, the inotropic response of cardiomyocytes to β-ARS is reduced [17,18], while the propensity to arrhythmogenic events increases. β-adrenergic response is also affected by ageing, and an age-dependent impairment between β-ARS and cardiac function has been demonstrated in both healthy and failing hearts [19,20]. β-ARS is altered as well in severe congenital heart disease patients [21]. Finally, recent studies also suggest a hyperactivation of the positive response of β-ARS patients with coronavirus disease 2019 (COVID-19), potentially leading to life-threatening arrhythmic events [22,23].
Refined knowledge of the role that each cellular component has within the β-ARS cascade, and of the consequences that may arise from disturbing its normal functioning, can therefore lead to a better understanding of different pathologies, as well as to the development of new pharmacological targets for their treatment. In these cases, mathematical modelling and simulation studies of β-ARS can be useful tools for investigating the mechanisms mediating arrhythmic events, assessing their multiscale consequences from the subcellular up to the organ levels, and designing effective treatments [24]. In this paper, we review the main roles of β-ARS in cellular cardiac electrophysiology, placing our emphasis on the description of the existing mathematical frameworks available for its representation and how insights obtained through experimental approaches have been integrated into these mathematical formulations. Then, we focus on how multiscale investigations exploiting mathematical descriptions of β-ARS have been used in different studies to better understand the role of β-ARS in heart disease and arrhythmias. We conclude our review by discussing a number of open questions requiring further combined experimental and computational investigations to improve our current understanding of the processes involved during β-ARS and their implications in cardiac function.

Mathematical Models of β-ARS
Several mathematical formulations, with varying degrees of complexity and physiological detail, have been proposed to date to describe different aspects of β-ARS in cardiac myocytes. In particular, special attention has been given to the modelling of β-ARS in mammalian ventricular electrophysiology, while studies on atrial electrophysiology are considerably lagging. Table 1 summarises the main β-ARS computational frameworks discussed in this section as landmark studies underlying the foundation of these modelling efforts. For the sake of simplicity, most models will be referred to by their first authors' names.
In general terms, cardiac myocytes present three different subtypes of β-adrenergic receptors (β-AR): β1-AR, β2-AR, and β3-AR [25], the former two being the most prevalent and important. Many of the existing mathematical models of β-ARS describe the bind-ing of an agonist, usually isoproterenol (ISO), to β1-AR and β2-AR and the subsequent phosphorylation cascade. In particular, the binding with β-ARs activates the receptorbound stimulatory G protein, which enhances the production of 3 -5 -cyclic adenosine monophosphate (cAMP), activating protein kinase A (PKA). PKA phosphorylates key cellular substrates that affect the functioning of several channels and ion pumps. The main targets located on the sarcolemma are the L-type calcium channels (I CaL ) [26], slow delayed rectifier potassium (I Ks ) channels [27], fast sodium current (I Na ) channels [28], the cystic fibrosis transmembrane conductance regulator (I CFTR ) [29,30], and the sodium-potassium pump (I NaK ) [31]. At the subcellular level, PKA also phosphorylates the ryanodine receptors (RyRs) and phospholamban (PLB) [26,32], both located on the sarcoplasmic reticulum, as well as troponin I (TPNI) [33], myosin binding protein-C [33], and titin [34], these located on the myofilaments. A graphical representation of the main PKA phosphorylation cascade is presented in Figure 1.  The first modelling approach for the consideration of β-ARS in cardiac electrophysiology described the effects of β-ARS by upscaling the magnitude of the most significantly upregulated ion channels during the β-adrenergic response (notably I CaL and I Ks ) or by shifting the activation curves of these currents [35,43]. Despite its simplicity, such an approach is sufficient to replicate to a good extent the main steady-state effects of β-ARS at the cellular level, such as action potential shortening, increased calcium transient amplitude, or a potentiated arrhythmogenicity. A more complex model developed by Greenstein et al. [37] accounted for these changes using a population-based approach, i.e., treating phosphorylation as a binary process independent of channel gating. This assumption implies that, at each time step, ionic concentrations and ionic currents or fluxes can be expressed as the sum of their phosphorylated and nonphosphorylated cellular parts. However, the steady-state nature of the abovementioned models imposes limitations for investigating transient behaviours in the activation of β-ARS or in dissecting the contribution to the β-ARS response of subcellular targets of PKA phosphorylation.
Saucerman et al. were the first authors to develop, in the early 2000s, a dynamic formulation of β-ARS [36,44]. The model incorporated a complete phosphorylation cascade of the L-type calcium channel and PLB in a rat ventricular myocyte electrophysiological model. The authors considered an extensive validation against published experimental data, including the temporal response of cAMP to ISO [45], PKA activation levels as a function of the concentration of cAMP [46], and PLB phosphorylation to ISO [47], together with experimental recordings of whole-cell patch-clamp I CaL current, calcium transients, and action potentials [48]. This seminal model of β-ARS has been expanded in multiple subsequent studies. Soltis and Saucerman added dynamic phosphorylation by Ca 2+ /calmodulindependent protein kinase II (CaMKII), in combination with the previously described PKA phosphorylation, to rabbit ventricular cardiomyocytes [29]. Phospholemman phosphorylation was also included in later works [49]. The model was also expanded by Negroni et al. to present an integrated framework of β-ARS and cardiac mechanical contraction [30]. More recently, Meyer et al. investigated a quasi-steady-state approximation of the Soltis and Saucerman model, demonstrating that their reduced-order formulation captured many of the features of the complete model while providing an efficient approximation over a broad range of parameter conditions [50].
The next landmark in β-ARS modelling was provided by Iancu et al., who proposed a new model considering cellular compartmentation of cAMP [38]. This led to a β-ARS model response dependent on local cAMP and PKA concentrations in the different cellular domains. The model response was compared to published measurements of cAMP activity under β-ARS [51]. The compartmentation introduced by the model provided a mechanistic explanation of how the high concentrations of cAMP as measured in experiments can modulate PKA activity.
The ideas described in the Saucerman and Iancu models were further developed by Heijman et al. [39]. In their work, they developed a compartmental model of β-ARS of the canine ventricular myocyte, including a dynamic β-ARS response integrated with CaMKII signalling. The main novelties of the model were the consideration of β1 and β2 isoforms in β-ARS stimulation and the use of a population-based approach with four different populations (phosphorylated by PKA, phosphorylated by CaMKII, phosphorylated by both, and nonphosphorylated) to represent the phosphorylation of its targets. The model was calibrated and validated using a wide range of canine and rodent experimental data, from measurements of molecular signalling (cAMP levels, adenylyl cyclase activity, PKA substrate [52,53]) to the complete cellular response (action potential duration, calcium transients [54]). The Heijman model has also been modified in later works, especially in terms of adjusting its underlying conductances and time constants to mimic the β-ARS response in other species, including humans [55,56].
A later model of β-ARS in mouse myocytes was developed by Bondarenko [40]. The model included some of the elements of the Heijman model, such as the compartmentalisation of the β-ARS system, but considered two types of L-type calcium channels and a different localisation of the RyRs in the compartmentalisation. The model matched experimental data from mice and showed how the different locations of some PKA phosphorylation targets can significantly affect action potential duration and calcium transients.
Recently, Khalilimeybodi et al. [41] formulated a model for β-adrenergic-induced cardiac hypertrophy in rat myocytes. The model included new molecular pathways that ended with the phosphorylation of the GSK3β and ERK molecules, which mediate gene expression in cardiac hypertrophy. The modelling results highlighted the significant contribution of these pathways in regulating cardiac hypertrophy in rats.

Differences in Formulation between Leading Mathematical Models of β-ARS
From the aforementioned models, the most widely adopted ones in modelling and simulation studies of β-ARS (either in their original form or with slight modifications) are the Soltis and Saucerman [29] and the Heijman [39] models. In these models, the biochemical reactions are modelled using ordinary differential equations or algebraic equations in cases with very fast reaction dynamics. The concentration of ISO is usually treated as an input to the model. In terms of complexity, the Soltis and Saucerman model is formed of 17 ordinary differential equations, containing 6 different PKA substrates (see Table 1). The Heijman β-ARS signalling model, on the contrary, consists of a coupled system of 57 ordinary differential equations and contains 8 different PKA substrates. The counting on these equations increases after the coupling of the signalling model with the respective cellular electrophysiology model, as illustrated below.
Both models share similar formulations of their phosphorylation targets, following a Michaelis-Menten formalism. To better illustrate differences in mathematical treatment, we consider the case of the inhibitory troponin subunit (TnI) as one of their PKA targets. Physiologically, the phosphorylation of TnI can increase the rate of calcium unbinding to calmodulin by 50% [57], Km trp being the binding constant. In the Soltis and Saucerman model, the time evolution of TnI phosphorylation is given by: where Km PKA and Km PP2A are Michaelis constants of phosphorylation and dephosphorylation, and K PKA and K PP2A are the corresponding phosphorylation and dephosphorylation rates of TnI. PKAC is the concentration of the PKA molecule, which is calculated by similar differential equations. TnI p and TnI np represent the concentrations of the phosphorylated and nonphosphorylated TnI, respectively. These are linked by conservation of species as: where the total TnI concentration, TnI total , is constant (70 µM). The phosphorylated PKA fraction is then given by the experimentally informed algebraic relationship: where f p TnI = TnI p /TnI total , and f po TnI = 0.0031. The phosphorylated fraction by PKA, f PKA TnI , can then be used to calculate the effective troponin affinity as: which will inform the respective differential equations to update the intracellular calcium balance.
The Heijman model also uses a Michaelis-Menten formalism for describing the phosphorylation rate, given by: which, instead of using concentrations, is written directly in this case in terms of phosphorylated and nonphosphorylated fractions, satisfying f p TnI + f np TnI = 1. The phosphorylated fraction by PKA is given in this model by the algebraic expression: with the effective troponin affinity finally calculated as: More specifically, Equation (7) exemplifies the population-based approach used in the Heijman model for each of its phosphorylation targets. The total contribution of a target is calculated as the weighted sum of two target populations (phosphorylated and nonphosphorylated). In Equation (7), these populations are described as Km np trpn and Km P trpn , each one with different kinetics (or conductances and time constants, in case of ion channels). Another important difference in the formulation between the two models can be seen in Equation (5). Following the compartmental nature of the model, the concentration of the free PKA catalytic subunit differs between the different considered compartments (caveolar, extracaveolar membrane, and cytosol), the concentrations in the latter (PKAC CYT and PP2A CYT ) being the ones affecting this particular process. Figure 2A illustrates the time evolution of the different PKA phosphorylation targets between two of the most recent versions of the two described models: the Xie et al. [58] version of the Soltis and Saucerman model and the Gong et al. [56] version of the Heijman model. The Xie version includes a new phosphorylation target, phospholemman (PLM), which affects the I NaK current, and a different time constant in I KS current phosphorylation. The Gong version is an adaptation of the Heijman model to human myocytes. Both models were integrated into the O'Hara and Rudy cardiomyocyte model [59] as one of the latest representations of human ventricular electrophysiology and paced at a basic cycle length of 1 Hz for 250 s, with activation of β-ARS by simulated ISO administration after 50 beats. The adaptation of the Xie model to human cardiomyocytes was conducted as in Pueyo et al. [60]. Considerable differences can be appreciated between the two models, the levels of target phosphorylation being higher in the Heijman model. Nevertheless, both models show similar trends in kinetics, with I KS presenting slower phosphorylation than I CaL or PLB and TnI being the ones with faster phosphorylation. The time evolution of the action potential duration exhibited similar behaviour in both models ( Figure 2B), with a small increase following ISO administration, followed by a reduction of approximately 50 ms after 250 s.
The differences between both β-ARS models at the level of ion channel currents and intracellular calcium fluxes are shown in Figure 3. As can be seen, the action potential shortening is similar in both models. However, the distinct formulations and phosphorylation levels can yield very different effects in specific currents and fluxes. This is for example the case for I CaL , I Nak , J rel , or J up , all exhibiting higher amplitudes when calculated with the Gong model, resulting in a considerably larger amplitude of its calcium transient.   [59] and after the simulation of the β-ARS response using the Gong [56] (yellow, dashed line) and the Xie [58] (blue, dotted lines) models. Results are presented for the last beat under the pacing protocol described in Figure 2. The ionic currents and intracellular calcium fluxes shown are: L-type calcium current (I CaL ); fast sodium current (I Na ); late sodium current (I NaL ); sodium-calcium exchanger (I NaCa ); sodium-potassium pump (I NaK ); rapid delayed rectifier potassium current (I Kr ); slow delayed rectifier potassium current (I Ks ); transient outward potassium current (I to ); calcium uptake via the SERCA pump (J up ) and ryanodine receptor calcium release (J rel ).

Multiscale Studies including β-ARS
The mathematical descriptions of β-ARS discussed before have been integrated into a variety of models of cardiac electrophysiology. Computational studies using the resulting multiscale cellular models have allowed quantitative investigations of the specific role of β-ARS in cardiac electrophysiology and different cardiac pathologies.
As a seminal example in the field, one of the grounds under the development of the earliest models of β-ARS was to resolve the ionic mechanisms underlying the formation of proarrhythmic triggers, such as EADs. The first model of β-ARS, developed by Zheng and Rudy [35], uncovered the nowadays well-established causative link between EADs and the reactivation of the L-type calcium current during the plateau and repolarisation phases of the ventricular action potential. Such an influential finding was made possible not only through mathematical model development but in combination with an exhaustive simulation study evaluating the effects of varying conductance amplitudes of different transmembrane currents, such as the L-type calcium channel or the potassium repolarising channels. A different approach to investigating the generation of EADs was presented by Tanskanen et al. [43], where they evaluated EADs formation using a β-ARS canine model with stochastic gating in the L-type calcium channels. Results suggested that the random fluctuations in channel gating can be the mechanisms underlying EAD generation.
The formulation of more complex models of β-ARS, such as those by Soltis and Saucerman [29] or Heijman et al. [39], allowed for a deeper physiological understanding of the underlying mechanisms. For example, Xie et al. [58] further studied the onset of EADs in rabbit ventricular myocytes using the Soltis and Saucerman β-ARS model, where each PKA target can exhibit different phosphorylation kinetics. As a result, they showed that the formation of EADs can be explained by the difference in phosphorylation levels between the L-type calcium and potassium channels. EAD generation in human myocytes under β-ARS was also studied in later works by Dai et al. [61]. Nevertheless, the range of potential applications in computational studies of β-ARS is certainly not limited to investigations of EADs as proarrhythmic triggers, as demonstrated below.

Effects of β-ARS in Electromechanical Coupling
Several computational studies have been conducted to evaluate the electrophysiological effects of β-ARS in the excitation-contraction response of cardiac myocytes. In particular, Kuzumoto et al. [62] studied the inotropic response using a guinea-pig ventricular model based on the Soltis and Saucerman model, demonstrating the positive inotropic effect of the I CFTR current under β-ARS. A later study by Negroni et al. [30] also developed a model of β-ARS that included contraction in rabbit myocytes. The model has also been used in [63] to evaluate the role of the PKA phosphorylation in additional myofilament targets and to study their impact on contractility, underscoring the importance of myofilament calcium sensitivity and cross-bridge cycling rate.
Recently, a study by Lyon et al. [64] presented a coupled model of human ventricular cellular electrophysiology [59] with a sarcomere contraction model [65] and investigated the effects of varying degrees of β-ARS in myocyte calcium and developed force during exercise. Among other findings, they identified synergic positive inotropic effects and differential lusitropic effects of β-ARS and sarcomere length on calcium and force dynamics.

Long QT Syndrome
Mutations in KCNQ1, the gene encoding the α-subunit of I Ks channels, are responsible for long QT syndrome type 1, the most common subtype of this congenital cardiac syndrome responsible for up to 30 to 35% of all cases. Most LQT1 patients experience cardiac events during episodes of β-ARS hyperactivity, such as physical exercise or acute emotional reactions [66]. The regulation of the I Ks channel activity by β-ARS has been studied in different works. Terrenoire et al. [67] studied it using the Soltis and Saucerman model, showing that alterations of kinetics in I Ks are critical in the control of action potential duration. Later studies also studied the role of β-ARS in long QT syndrome using the Heijman model in human myocytes [55,68]. O'Hara and Rudy [55] showed that under the mutation Q357R of the KCNQ1 gene, arrhythmias are caused by the reduced I Ks channel expression, whereas Kang et al. [68] found that β-ARS increased I KS contribution to repolarisation. This can lead to an augmentation of the transmural dispersion of repolarisation, which can be arrhythmogenic.

Heart Failure
Heart failure patients are characterised by a long-term remodelling of multiple cellular processes, including calcium handling, excitation-contraction coupling, and betaadrenergic signalling [69]. These alterations can lead to arrhythmia, as studied by Morotti et al. [70] through the development of a mouse ventricular model integrating β-ARS based on the Soltis and Saucerman model. The model simulations showed that during heart failure conditions, an increase in sodium concentration can generate a positive feedback loop with CaMKII, increasing phosphorylation of CaMKII targets, thus raising intracellular calcium. Such a positive feedback loop is highly arrhythmogenic, and the study demonstrated it can be potentiated by the β-ARS response. Gong et al. [56] recently adapted the Heijman model to human ventricular electrophysiology and systematically explored what alterations in the cell targets can give place to the phenotypic characteristics of β-ARS in the disease. One additional strength of their work was the consideration of the variability in model parameters, reflecting the heterogeneity between different experimental measurements or individuals, which can greatly affect the simulation outcomes in terms of proarrhythmic findings [71,72]. The model by Gong et al. was also recently used in a study by Mora et al. [73], where they reproduced the electrophysiological response to β-ARS in heart failure and analysed changes in the electrophysiology response. They found that the inotropic effect of β-ARS was reduced in heart failure due to altered calcium handling and that these myocytes have a higher probability of developing EADs. Their simulations also showed that the response to β2-ARs was blunted in heart failure, with β1-ARs being the biggest regulators of most electrophysiological effects under β-ARS.
A different approach in the modelling of heart failure is the one presented by Loucks et al. [74]. In their work, they included the detubulation process (loss of T-tubules of the cell membrane) that is usually concomitant to the remodelling of cardiac myocytes in heart failure. This process affects the distribution of the L-type calcium channels, a target of the β-ARS phosphorylation cascade. The effects of the combination of the detubulation and β-ARS in heart failure and its arrhythmogenicity were also studied in a later study by Sanchez-Alonso et al. [75], showing that redistribution of the L-type calcium channels to the sarcolemma can alter myocyte function and trigger pathological events.

Other Cardiac Pathologies
Mathematical models of β-ARS have been also used in combined experimental and computational studies to investigate the generation of DADs in catecholaminergic polymorphic ventricular tachycardia [76], the suppression of alternans in the postinfarction border zone [77,78], and the occurrence of re-entrant arrhythmias in Brugada syndrome [79] or to study the possible proarrhythmic effects of two anti-COVID-19 drugs under β-ARS [23]. Finally, the effects that low-frequency oscillatory patterns (e.g., respiration) can have in the sympathetic nervous system, and therefore in the β-ARS response of human heart function, have been experimentally characterised and computationally simulated in several recent works [60,80].

Discussion and Perspectives
The importance of β-ARS in regulating cardiac function, as well as its role in promoting arrhythmic events, is nowadays well established. During the last 20 years, a substantial number of modelling and simulation studies have greatly contributed to consolidating the present understanding of the cellular mechanisms underlying the multifactorial response to β-ARS [5,16,20,81,82]. In fact, the changes that β-ARS mobilises at different heart scales can considerably affect organ performance, and the different studies covered in this review paper have underlined the need of including β-ARS to advance the characterisation of (patho-)physiology. It is for these reasons that further advances on β-ARS modelling, and its wider integration in simulation studies of heart function, are expected to considerably advance the potential of these technologies towards precision medicine and the understanding of diseased states, especially in investigations of arrhythmogenicity, exercise, or mechanical response.

Future Work in β-ARS Modelling
One of the main challenges ahead in the modelling of the β-ARS response and function is to fully understand cAMP and PKA signalling compartmentation. Physiologically, cellular compartmentation is produced by the contribution of different processes that concentrate some molecular transit in confined spaces within the cell [83]. In that way, specific cellular functions can be performed more efficiently. It has been shown that different pathologies such as hypertrophy or heart failure can alter cell compartmentation, affecting cardiac function [84]. In combination with new experimental technologies such as fluorescence resonance energy transfer (FRET) [45,63], additional modelling efforts remain necessary to reliably reflect altered compartmentation in mathematical electrophysiological models and to uncover its functional consequences on cAMP and PKA signalling.
Further experimental characterisation of phosphorylation time dynamics is also expected to notably contribute to the validation and refinement of mathematical models of β-ARS. Both experimental [85] and modelling studies [58,86] have demonstrated that differences in the time course of activation under ISO between different PKA targets can lead to proarrhythmic triggers in the form of EADs. Beyond the studies on compartmentation suggested above, technologies such as FRET also hold great potential to further inform the mathematical modelling of β-ARS in terms of dynamic changes in phosphorylation levels, whose functional consequences are still not clearly understood [87]. Moreover, the extension of β-ARS models to incorporate additional phosphorylation targets, such as titin and myosin binding protein-C, would also qualify such models for quantitative investigations on how their active role in cardiac contraction and electrophysiology is modulated by sympathetic stimulation.
Some of the mechanisms behind how β-ARS affects the normal functioning of ion channels remain unclear. For example, during chronic β-ARS, I Ks presents the opposite effect (reduction of I Ks current) that under acute β-ARS (I Ks increase) [88]. Conditions of chronic β-ARS also pertain to different diseased states and are for instance characteristic of heart failure. Furthermore, several studies [89,90] have suggested a possible mismatch between calcium-dependant and voltage-dependant inactivation of the I CaL channel under β-ARS, potentially leading to arrhythmia formation.
The scarcity of human electrophysiological data implies that the existing mathematical models of β-ARS are largely based on data from other species then carefully reparametrised to reproduce the behaviour observed in human cardiomyocytes [55,56]. Apart from the well-known differences in ion channel conformation and distribution [91], differences have also been reported in the molecular processes involved in β-ARS and in the overall response of human cardiomyocytes to adrenergic stimulation, presenting opposite effects in the shortening of the action potential duration or in the functioning of the sodium-calcium exchanger [90,92,93]. In addition, the density and anatomical location of β-ARS receptors in the heart are different between species [53,94]. Further data acquisition and model developments will allow ascertaining the implications of these interspecies differences in β-ARS in heart function.

β-ARS and 3D Modelling
The β-ARS response is present at the whole-organ level, including both the ventricles and the atria. For example, many studies have linked β-ARS with the occurrence of episodes of atrial fibrillation [5,95], and in fact, beta-blockers are first-line therapy for rate control in atrial fibrillation patients [96]. However, most modelling efforts of β-ARS to date have been placed on ventricular electrophysiology, while only a small number of recent studies have addressed its modelling in atrial cells, specifically in mice [97,98]. Extending the existing modelling frameworks of β-ARS to human atrial electrophysiology, from the cellular up to the 3D whole-atria level, could lead to significant advances in understanding the involvement of the sympathetic nervous system in atrial arrhythmias, as done in the past in terms of parasympathetic stimulation [99][100][101].
Moreover, β-adrenergic receptors are heterogeneously distributed along different areas of the heart, especially under pathological conditions [81]. Transmural and apicobasal heterogeneities in the sympathetic innervation system have also been described [102,103]. In addition, through its regulation of fast sodium channels, β-ARS also affects the conduction velocity of myocardial tissue [68,104]. The degree to which these spatially distributed regions of heterogeneous action potential and conduction velocity can contribute to providing suitable arrhythmogenic substrates can only be fully realised by upscaling the modelling and simulation efforts of β-ARS to multiscale whole-ventricular, whole-atria, or whole-organ levels. In particular, the resources provided by the "Stimulating Peripheral Activity to Relieve Conditions" (SPARC) programme, including tools and curated and annotated high-resolution datasets in different species for investigations on autonomic nervous stimulation (Figure 4), may prove fundamental in understanding the role of a heterogeneous innervation in heart function [105]. However, there is a lack of computational studies to date implementing β-ARS models in 3D geometries. One of the reasons is the increased complexity of these models, which escalate the size of the state space by splitting the β-adrenergic targets into phosphorylated and unphosphorylated counterparts. Nevertheless, the main reasons underpinning their reduced uptake in 3D studies might actually come from the enormous computational cost needed to conduct simulations in time scales of relevance for β-ARS effects, usually after several tens to hundreds of heartbeats. Current whole-heart human 3D models are often conformed by meshes of millions of elements, and a vast computational power would be required to solve all model equations on them for such long time scales [24,106].
The use of advanced mathematical techniques can play a significant role in the upscaling of β-ARS models to the whole-organ level. The recent work by Meyer et al., applying model reduction techniques by quasi-steady-state approximations, was able to simplify the β-ARS cascade of the Soltis and Saucerman model into a minimal version of only two ordinary differential equations while retaining most of the relevant features of the original model [50]. Another approach, already used in the study by Sanchez-Alonso et al. [75], is the direct implementation of the β-ARS effects in the model. The changes in the magnitude of the ionic currents that were observed in the Heijman model after PKA phosphorylation were applied to the electrophysiological model used for the 3D simulations. The computational cost is lower, as there is no coupling with the β-ARS model, and thus the number of equations that need to be solved is reduced. However, such an approach is suboptimal, as part of the changes induced by β-ARS (such as shifts in activation/inactivation curves, or the influence on gate activation kinetics) is effectively lost. A better solution could be the full implementation of the coupled model (electrophysiology + β-ARS), with a 3D model initialisation based on the state variables of the cellular model after the appropriate number of beats (depending on the level of β-adrenergic phosphorylation of interest). This initialisation approach, widely exploited in other 3D electrophysiological and electromechanical studies [107,108], could allow the simulation of different stages of the β-ARS cascade in whole 3D models.

Conclusions
Over the last two decades, great advances have been made towards the mechanistic modelling and multiscale simulation of the cardiac β-ARS response. While significant progress has already been accomplished, more complete and comprehensive descriptions of this physiologically crucial signalling process are expected to further contribute to an improved understanding of its role in healthy and pathological conditions. In particular, the integration of these β-ARS models into whole-heart computational models, linked with tissue electrophysiology and mechanics, is envisioned as the next crucial step in the realisation of their full potential to bring novel insights into the modulation of human cardiac function. These perspectives will qualify digital twin technologies with enhanced capabilities for arrhythmic risk stratification under physical exercise or emotional stress, as well as for the safety pharmacology and personalisation of therapies, such as beta-blockers, one of the most recurred treatments in patients of heart failure, atrial fibrillation, survivors of myocardial infarction, or at risk of suffering major adverse cardiac events.  Data Availability Statement: The codes used in this work can be found at: https://github.com/ JQXGong/Ohara-beta-adrenergic (Gong et al. [56]) and https://github.com/saucermanlab/modelarchive (Soltis and Saucerman [29]) (accessed on 4 July 2021).