A multiscale predictive digital twin for neurocardiac modulation

Cardiac function is tightly regulated by the autonomic nervous system (ANS). Activation of the sympathetic nervous system increases cardiac output by increasing heart rate and stroke volume, while parasympathetic nerve stimulation instantly slows heart rate. Importantly, imbalance in autonomic control of the heart has been implicated in the development of arrhythmias and heart failure. Understanding of the mechanisms and effects of autonomic stimulation is a major challenge because synapses in different regions of the heart result in multiple changes to heart function. For example, nerve synapses on the sinoatrial node (SAN) impact pacemaking, while synapses on contractile cells alter contraction and arrhythmia vulnerability. Here, we present a multiscale neurocardiac modelling and simulator tool that predicts the effect of efferent stimulation of the sympathetic and parasympathetic branches of the ANS on the cardiac SAN and ventricular myocardium. The model includes a layered representation of the ANS and reproduces firing properties measured experimentally. Model parameters are derived from experiments and atomistic simulations. The model is a first prototype of a digital twin that is applied to make predictions across all system scales, from subcellular signalling to pacemaker frequency to tissue level responses. We predict conditions under which autonomic imbalance induces proarrhythmia and can be modified to prevent or inhibit arrhythmia. In summary, the multiscale model constitutes a predictive digital twin framework to test and guide high‐throughput prediction of novel neuromodulatory therapy.

r A key feature of the neurocardiac computational model is the connection between the autonomic nervous system and both pacemaker and contractile cells, where modification to pacemaker frequency drives initiation of electrical signals in the contractile cells.
r We utilized atomic-scale molecular dynamics simulations to predict the association and dissociation rates of noradrenaline with the β-adrenergic receptor.
r Multiscale predictions demonstrate how autonomic imbalance may increase proclivity to arrhythmias or be used to terminate arrhythmias. r The model serves as a first step towards a digital twin for predicting neuromodulation to prevent or reduce disease.

Introduction
Every affection of the mind that is attended with either pain or pleasure, hope or fear, is the cause of an agitation whose influence extends to the heart. William Harvey (1628) Modulation of peripheral nerve activity holds great promise for treating, preventing and reversing a host of cardiovascular diseases, including heart failure, hypertension and heart rhythm disorders (Fukuda et al., 2015;Grassi, 2007;Grassi et al., 1998Grassi et al., , 2007Grassi, Seravalle et al., 2008;Lahiri et al., 2008;Lopshire et al., 2009;Randall & Ardell, 1985;Zuanetti et al., 1987). Even though the promise of this approach is well accepted, successful design and implementation of neuromodulatory therapy is limited by the utter paucity of tools that can predict cardiac and vascular responses to peripheral nerve activity. To develop such tools, a comprehensive collaborative effort to map and model complex spatio-temporal subcellular, cellular and tissue-level responses is needed. Here, we present the first neurocardiac simulator resulting from a collaborative team as part of the NIH Common Fund program Stimulating Peripheral Activity to Relieve Conditions (SPARC).
The computational model we present was built using multiscale data to both inform and validate the model parameters and outputs both within subcellular, cellular and tissue scales as well as across them. To our knowledge, this constitutes the first instance of a digital twin for neurocardiac modulation. The tool can be used to reveal how the autonomic nervous system regulates cardiac function and to predict neuromodulation to prevent or reduce disease. The simulator integrates data from the atomic (e.g. β-adrenergic receptor (βAR)-noradrenaline (NA) interaction affinities and rates), subcellular (e.g. cAMP dynamics), cellular (e.g. action potential and calcium transient), tissue (e.g. conduction velocity) and system (e.g. pseudo-ECG) scales. The model predicts the effect of efferent stimulation of the sympathetic and parasympathetic branches of the autonomic nervous system (ANS) on the cardiac sinoatrial node (SAN) and ventricular myocardium and accounts for the coupling between them.
Development and application of computational modelling and simulation in the investigation of plausible fundamental cardiac mechanisms and to predict cardiac physiology and pathophysiology has grown dramatically over the last two decades. Indeed, multiple modelling studies have begun to address varying aspects of autonomic cardiac control. As experimental and clinical approaches have become increasingly sophisticated, new data have revealed multiple physiological processes and mechanisms from the protein to the patient. Computational modelling and simulation have kept pace with these developments and a number of important models are now available in the literature and in the public domain.
Computational and mathematical models have also been used to demonstrate functional mechanisms of heart rate variability (Karavaev et al., 2019;Prokhorov et al., 2021;von Rosenberg et al., 2019), applied to reveal the dynamics of autonomic control during orthostatic adaptation (Ishbulatov et al., 2020) and shown to capture heart rate control by sympathetic and parasympathetic discharge (Kember et al., 2011(Kember et al., , 2017Prokhorov et al., 2021;Warner & Cox, 1962). Detailed subcellular models have been utilized to predict the mechanisms of postsynaptic receptor activation following nerve stimulation and to propose likely mechanisms for feedback and feed forward regulation of cardiac electrophysiology and contraction in the SAN and in contractile cells, respectively (Behar et al., 2016;Castellanos & Godinez, 2015;Iancu et al., 2007Iancu et al., , 2008O'Hara & Rudy, 2012;Yang & Saucerman, 2012).
There have also been numerous atomistic modelling and simulation studies of βARs that link the autonomic nervous system to downstream cellular signalling events and function recently reviewed (e.g. Hilger, 2021;Ribeiro & Filizola, 2019;Wang & Miao, 2019). These studies were enabled in large part by advances in atomic-resolution structural characterization of adrenergic receptors and their complexes with small-molecule ligands as well as regulatory proteins via X-ray crystallography, cryogenic electron microscopy, NMR spectroscopy and other experimental techniques (Katritch et al., 2013;Kobilka, 2011;Shimada et al., 2021;Thal et al., 2018). Exploring atomic-resolution receptor and ligand conformational dynamics was found to be crucial to provide an accurate description of βAR signalling pathways (Hilger, 2021) but accurately determining ensuing energetic and kinetic quantities was proven to be challenging and requires usage of simplified models, extremely long simulations and/or advanced simulation techniques. This is, however, essential for elucidation of atomistic simulation-derived parameters for functional kinetic modelling of subcellular signalling such as those from the Iancu-Harvey (Iancu et al., 2007) and Soltis-Saucerman (Soltis & Saucerman, 2010) models as used in our study.
The frequency of activation of pacemaking cells in the SAN, the conduction velocity of action potentials, as well as the strength of excitation-contraction coupling, is tightly regulated by the autonomic nervous system. However, despite the importance of the fundamental connection between the brain and heart, there is still an incomplete understanding of how the autonomic nervous system regulates cardiac function in health and disease (Armour, 2004;Olshansky, 2005;Patterson et al., 2007;Tan et al., 2008;Verrier & Antzelevitch, 2004;Young et al., 1993;Zhou, Jung et al., 2008;. Description and prediction of the mechanisms underlying the interaction between nervous system discharge and the resultant emergent cardiac events in a digital twin are needed to allow for identification and specific targeting to prevent and treat arrhythmia-provoking conditions by drugs and by direct electrical stimulation (Hanna et al., 2021;Meng et al., 2017).

Autonomic nervous system network structure
The structure of the ANS network model is based on the connectivity among various intrathoracic neurons and central neurons that are involved in cardiac regulation as described by Shivkumar et al. (Fukuda et al., 2015;Shivkumar et al., 2016) and the corresponding computational models of Kember et al. (Kember et al., 2011) (see Fig. 1). The neural circuit consists of the sympathetic pathway (SNS) and the parasympathetic pathway (PNS). The SNS comprises three layers (subnetworks): the central nervous system (CNS), the intrathoracic nervous system (ITNS) and the intrinsic cardiac nervous system (ICNS). Efferent pathways connect the CNS to the ITNS, and the ITNS to a subpopulation of neurons in the sympathetic ICNS (i.e. the S-ICNS). Afferent pathways connect the S-ICNS with the ITNS. The modular workflow consists of a network layer representation for the sympathetic branch (blue) and the parasympathetic branch (green) of the autonomic nervous system (ANS) that both synapse onto the cardiac sinoatrial node (SAN) as well as on cardiac ventricular myocardium (purple box). The model allows for prediction of efferent sympathetic and parasympathetic signalling and target-organ responses. The model drives spontaneous pacemaker action potentials modulated by the ANS that determine the cardiac ventricular pacing frequency (indicated by the yellow arrow). [Colour figure can be viewed at wileyonlinelibrary.com] The PNS comprises two layers (subnetworks): the CNS and a subpopulation of neurons in the parasympathetic intrinsic cardiac nervous system (i.e. the P-ICNS), with efferent connections from the CNS to the P-ICNS. The S-ICNS and the P-ICNS are also interconnected. Note that there is no clear delineation between a sympathetic ICNS and a parasympathetic ICNS, and the S-ICNS and the P-ICNS in the model are extensively connected, effectively forming a single subnetwork.

ITNS and ICNS neuronal dynamics
The ITNS and ICNS subnetworks each contain 50 model neurons. Individual neurons in the ITNS and ICNS are described by a single-compartment conductance-based generalized integrate-and-fire (GIF) model (Jolivet et al., 2004). The subthreshold dynamics of the neurons include a leakage current, a delayed-rectifier potassium current, an M-current (Kwag et al., 2014) and a synaptic current. The equations governing the subthreshold dynamics of the j-th neuron are: where t is time (in ms), v j is the transmembrane potential (in mV); x j = n j , w j are the activation variables of the delayed rectifier potassium current and the M-current, respectively; I syn, j is the input synaptic current into the j-th neuron; C m is the membrane capacitance, and g L , g K , g M are the maximal conductances of the leakage, delayed-rectifier potassium and M-currents, respectively; E L , E K , E M are the reversal potentials of the leakage and delayed-rectifier potassium and M-currents; τ k = τ n , τ w are the time-constants for the activation variables, and the function x ∞ (v) is the subthreshold steady-state activation of the activation variable x j , given by the equation: When the transmembrane potential of a model neuron reaches a threshold potential (v T ), the neuron fires an action potential and the transmembrane potential is 'reset' . Note that in the GIF model action potentials are not explicitly modelled. Each action potential can be thought of as an instantaneous spike in membrane potential followed by a brief absolute refractory period (t ref ). Following the absolute refractory period, the transmembrane potential is set to a reset potential (v reset ), and the subthreshold dynamics (described above) are resumed. To capture the effect of the action potential on the gating variables, n and w are updated to include the increase that would occur during a stereotypical spike. Specifically, when v j (t ) = v T , the membrane potential and gating variables are updated to: To model the saturation of the gating variables (i.e. n and w must be between 0 and 1), the update values for the gating variables, n and w, are scaled by (1 − x j (t )): . GIF model parameters where set to elicit similar firing properties of autonomic neurons (Springer et al., 2015) (see Fig. 2) and are provided in Table 1. Except for the M-current conductance g M and synaptic conductances, parameters for the neurons were homogeneous; g M was chosen from a beta distribution with minimum and maximum values taken to approximate values from previous models (Kwag et al., 2014). Synaptic dynamics are described below.

CNS neuronal dynamics
The sympathetic and parasympathetic CNS subnetworks contain 50 model neurons each. Spike times for the individual neurons in the CNS subnetworks are modelled as renewal processes. That is, the firing times of the j-th neuron in the sympathetic and the parasympathetic branches of the CNS are given by Poisson processes with firing rates λ SNS and λ PNS and a post-spike absolute refractory period of t ref (Gerstner & Kistler, 2002). In simulations when the SNS is activated, λ SNS = 0.5 Hz, and when the PNS is activated, λ PNS = 0.25 Hz; otherwise they are 0.

Synaptic dynamics
The synaptic dynamics throughout the neural circuit are modelled as alpha function synapses (Gerstner & Kistler, 2002;Lewis & Rinzel, 2003). It is assumed that each time a presynaptic neuron fires there is a stereotypical increase and decrease in the synaptic conductance of the post-synaptic neuron of the form: whereḡ syn sets the maximal conductance, τ r is the rise time constant and τ d is the decay time constant. When the presynaptic neuron fires multiple times the stereotypical increases and deceases in the post-synaptic conductance  (Springer et al., 2015) and (B) model predicted neuronal firing recorded by using a 1 s stimulus with differing stimulus current input amplitude. A, experimental recordings indicate that neurons generate at least three types of firing patterns: tonic, accommodating and phasic. Tonic neurons fired repetitively at frequencies proportional to strength of stimulus. Accommodating neurons adapted and ceased firing at lower stimulus levels. Phasic neurons fired one to four spikes and ceased firing. B, model neurons displayed tonic, accommodating and phasic firing dynamics. Variations in maximal M-current conductance (g M ) elicited a change in firing dynamics observed in experiments and indicate a plausible mechanism to explain the firing dynamics in A. (i.e. the damped harmonic oscillator). Therefore, because the synaptic responses are assumed to add linearly, when the j-th neuron fires at times t spike,k , its synaptic output g syn s j is given by: with update conditions: Values of the synaptic time constants (τ r andτ d ) were set to values estimated in Selyanko et al. (1979) and Wheeler et al. (2004).

Network connectivity
Within the ITNS and ICNS subnetworks of the SNS and within the ICNS subnetworks of the PNS, random recurrent connections from a presynaptic neuron j to a postsynaptic neuron k occur with probability p ITNS and p ICNS . Neurons are excitatory with probability p and inhibitory with probability 1 − p. The probability of connection between the S-ICNS and P-ICNS neurons is p SP . In the PNS, efferent connections from CNS neurons to ICNS neurons occur randomly with probability of connection p PNS . In the SNS, afferent connections from the ICNS to the ITNS occur with probability p SNS,aff . In the SNS, neurons in the ITNS receive efferent input from CNS neurons with probability p ITNS,eff , and ICNS neurons receive efferent input from ITNS neurons with probability p ICNS,eff . Each neuron in the ITNS that receives efferent input receives it from n+1 random presynaptic neurons in each preceding layer (with n = 3), where the n connections are weak, and the 1 connection is strong. This n+1 convergent innervation has been observed in both amphibian sympathetic ganglia (Wheeler et al., 2004) and mammalian sympathetic ganglia (Hirst & McLachlan, 1986;Janig & McLachlan, 1992;Skok & Ivanov, 1983). The random connectivity in the neural circuit inherently generates populations of afferent, efferent and local circuit neurons (Fukuda et al., 2015;Shivkumar et al., 2016).
To model afferent feedback from the cardiovascular system, each neuron in the ICNS and ITNS receives independent periodically modulated stochastic input. The input is described by a non-stationary Poisson process with arrival function μ(t ) that determines the firing times of the presynaptic neurons, synapsing onto the ITCNS and ICNS neurons. μ(t ) is a piecewise-linear function with a period of 1000 ms that mimics the shape of the arterial pressure during the cardiac cycle: , with time t in ms. The amplitude A of the input to ICNS neurons and ITNS neurons is 0.003 Hz and 0.002 Hz, respectively, so that the subnetworks that are anatomically closer to the heart receive stronger cardiac input (Kember et al., 2011). Note that feedback is applied only in simulations described in Fig. 3. The synaptic output of the cells is included in the synaptic current of each individual neuron, described by: where w k k, j and w C k, j are synaptic weights, and E syn, f (k) and E syn,C are the reversal potentials for the synaptic currents between neurons and the network input, respectively. The function f(k) outputs e or i, depending on whether the k-th cell is excitatory or inhibitory. The synaptic weights were drawn randomly from a distribution tuned using the model neurons without an M-current, so the fraction of phasic neurons in the neuronal network approximates those in the experimental data in Beaumont et al. (2013). These synaptic weights were also used to implement the n+1 convergent efferent innervation protocol in the SNS. The reversal potentials E syn,e and E syn,i were chosen to approximate the values represented by nicotinic receptor and GABA receptor activation, respectively.

Figure 3. Predicted distribution of firing types in the ANS layered network model
A, schematic representation of the autonomic nervous system layered network model. Each rectangle (CNS, ITNS, ICNS) represents a random network of model-generated integrate-and-fire neurons with delayed rectifier potassium, leak, M-type potassium and synaptic currents. Synaptic currents generate the intra-and inter-network connections. Inter-network connections are represented by arrows. Temporal correlation to left ventricular pressure is used to classify firing types as follower, phasic or tonic. Distribution of firing types within each network was predicted at steady-state, averaged over seven 60-s simulations. Model-generated neuronal firing patters that did not clearly fit into one classification were labelled as 'other'. B, examples of neural firing relative to cardiac phase in ICNS neurons measured in dog (Beaumont et al., 2013). Probability density of ICNS neuronal firing (i.e. ICNS firing histograms) as a function of the timing within the left ventricular pressure (LVP) cycle (represented as a phase between 0 and 2π ). The LVP is indicated by the grey curve. [Colour figure can be viewed at wileyonlinelibrary.com]

Functional models of cardiac cells and tissues
Sinoatrial node cell model representations. The layered network model of the autonomic nervous system generates both sympathetic and parasympathetic nerve stimulation outputs, which are used as scaling factor inputs for β-adrenergic stimulation and muscarinic stimulation, respectively in the Behar-Yaniv rabbit SAN model (Behar et al., 2016) to produce spontaneous action potentials that control rates of cardiac pacemaking. Ventricular cell model representations. We then merged the Iancu-Harvey model of cyclic adenosine monophosphate (cAMP) compartmentation (Iancu et al., 2007(Iancu et al., , 2008 with the Soltis-Saucerman ventricular cell electrophysiology model (Soltis & Saucerman, 2010), which includes all the relevant components required for a detailed analysis, including accurate cellular electrophysiology, Ca 2+ handling (Shannon et al., 2004) and the cAMP-dependent protein kinase (PKA) (Saucerman et al., 2003) and Ca 2+ /calmodulin-dependent protein kinase (CaMKII) phosphorylation pathways (Soltis & Saucerman, 2010). Pacing frequency was set to the heart rate generated by the Behar-Yaniv SAN model.

One-dimensional cable
One-dimensional (1D) ventricular heart tissue was simulated as a fibre of 165 cardiomyocytes (1.65 cm) (Glukhov et al., 2010) with reflective boundary conditions. Transmural heterogeneity was incorporated into the tissue by a linear decrease to maximal I Kr conductance (Myles et al., 2010), corresponding to a linear transition from epicardial to endocardial tissue (Fedida & Giles, 1991). The diffusion coefficient D x was set to 0.002 cm 2 /ms to establish a conduction velocity of 61-73 cm/s (epicardium-endocardium in wild-type (WT) conditions) (Brugada et al., 1990).

Noradrenaline binding affinity to the β-adrenergic receptor
Structural modelling of type 2 β-adrenergic receptor (β 2 AR) was performed with ROSETTA software (Bender et al., 2016;Yarov-Yarovoy et al., 2006) using its active-state crystal structure (PDB ID: 4LDO) as a template. The intracellular loop 3 (ICL3) region missing in the template structure was replaced by a three-residue linker predicted using Rosetta cyclic coordinate descent (CCD) loop modelling application (Wang et al., 2007). Ten thousand decoy structures were generated with the top one selected by total score for further simulations. ROSETTALIGAND (Davis & Baker, 2009;Meiler & Baker, 2006) was used for β 2 AR J Physiol 601.17 cationic (+) and neutral (0) noradrenaline (NA) docking calculations. Fifty thousand structures were generated with the top 10% selected by total score, out of which the one with the lowest interfacial score (interaction energy) was chosen. CHARMM-GUI (Jo et al., 2008) was used to prepare β 2 AR-NA complexes embedded in the palmitoyl-oleoyl-phosphatidylcholine (POPC) lipid bilayer solvated by a 150 mM NaCl aqueous solution. We used CHARMM36m (Huang et al., 2017) protein and C36 lipid (Klauda et al., 2010) force fields, TIP3P water model (Jorgensen et al., 1983) and NA(+) and NA(0) parameters based on CHARMM general force field (CGENFF) (Vanommeslaeghe et al., 2010). The latter were developed using the CGNEFF program  for initial parameter estimates, which were then optimized through their subsequent fitting to quantum mechanical (QM) reference data using an established protocol (Vanommeslaeghe et al., 2010). All-atom molecular dynamics (MD) simulations were run with NAMD (Phillips et al., 2005) in the NPT (constant number of particles, N, pressure, P, and temperature, T) ensemble at 310 K and 1 atm pressure. Initial equilibration MD simulations were run for ∼52 ns with the first ∼2 ns using gradually reduced harmonic restraints on the protein and lipid tail carbon atoms followed by unrestrained 50 ns-long MD simulation runs. They were used to start 200 ns long well-tempered metadynamics (Barducci et al., 2008) MD simulations to estimate ligand binding affinities. Distance between centres of mass (COM) of protein α-helical core and NA(+) was used as reaction coordinate R extending to 50Å in bulk, whereas angle θ between vectors connecting these points and ligand orthosteric binding site was used to restrain sampling within θ ≤ 30°cone as was done in Provasi et al. (2009) andSchneider et al. (2015). Affinities in the form of dissociation constants, K d , were computed as was done previously (Provasi et al., 2009;Schneider et al., 2015) using equation where is the solid angle defined by the conical restraint, R bulk = 50 Å, N A is the Avogadro number, W(R) is free energy profile, k is the Boltzmann constant and T is absolute temperature. NA K d value at pH 7.4 was computed from corresponding estimates for NA(+) and NA(0) using the Henderson-Hasselbach equation, as was done previously (Yang et al., 2020). NA(+) dissociation or 'off ' rate constant kr LR was estimated using Kramer's rate formalism using equation, (Dorairaj & Allen, 2007). The free energy profile from metadynamics MD simulation of the β 2 AR-NA(+) system was used to estimate activation free energy barrier height ( W activation ) and corresponding curvatures at free energy maxima (W (R barrier )) and minima (W (R well )).
The diffusion coefficient at the free energy maximum, D(R barrier ), was estimated using position autocorrelation function from last 20 ns of separate 30 ns long MD simulations with harmonic restraints on the NA(+) position as was done previously (Vorobyov et al., 2014). The z-component of a distance between COM of the protein α-helical core and NA(+) was used as reaction coordinate for diffusion coefficient calculations. NA(+) association or 'on' rate constant, kf LR , was then computed as kr LR /K d .
In our one-dimensional tissue model, the β-AR binding to ligand equation (3) (Yang & Saucerman, 2012) was used in the Iancu-Harvey model of cAMP compartmentation (Iancu et al., 2007;Iancu et al., 2008) for calculating concentration of ternary complex (LRG), ligand-receptor (LR) complex and receptor bound G-protein (RG). Gs is free G-protein concentration, and [Stimulus] is input stimulus from SNS and PNS. The rate constants k f LR and kr LR were obtained from molecular dynamics (MD) simulations as described above. dLR

Current density changes in diseased heart
To simulate diseased heart, we modified the current density changes based on heart failure conditions as shown in Table 3. In addition, CaMKII expression is increased in failing human myocardium (Hoch et al., 1999). We simulated CaMKII overexpression (CaMKII-OE) as in Soltis-Saucerman (Soltis & Saucerman, 2010). In all simulations, we used a second-order Runga-Kutta method with a time step of 1/128 ms (∼7.5 μs) and included a linear interpolation scheme to compute the spike times between the time steps in the generalized integrate-and-fire (GIF) model (Shelley & Tao, 2001). The interpolation scheme ensures that the method is truly second order, despite the discontinuities due to the fire-and-reset conditions in the GIF formulation. Code for simulations and analysis was written in C++ and MATLAB 2018a (The MathWorks, Natick, MA, USA). Code was run on an Apple Mac Pro machine with 2 2.7 GHz 12-Core Intel Xeon processors, and an HP ProLiant DL585 G7 server with a 2.7 GHz 28-core AMD Opteron processor. Code was compiled with the Intel ICC compiler, version 18.0.3. Numerical results were visualized using MATLAB R2018a. The models used in this paper are available on SPACR Portal (https://sparc.science/datasets/ 322?type=simulation).

Results
In this study, we endeavoured to develop a neurocardiac simulator suite containing modelling tools that will allow for prediction of a range of autonomic effects on cardiovascular function. The model framework is shown schematically in Fig. 1 as a modular workflow comprising a network layer representation of the sympathetic branch (blue) and the parasympathetic branch (green) of the ANS that both synapse onto the cardiac SAN as well as on cardiac ventricular myocardium (purple box). The model allows for the prediction of efferent sympathetic and parasympathetic signalling and model representations of target-organ responses.
We first set out to develop computational models of single neuron action potentials from the three-layer representation of the autonomic nervous system. In Fig. 2A, the firing dynamics from experimental recordings of neuronal action potentials measured from the superior cervical ganglion of adult male rats are shown (Springer et al., 2015). An increasing amplitude of current injection applied for 1 s is shown from the top to the bottom panels. Neurons displayed three distinct notable types of firing: tonic, accommodating and phasic. Tonic neurons fired repetitively at frequencies related to strength of stimulus. Accommodating neurons adapted and demonstrated reduced firing at lower stimulus amplitudes. Phasic neurons fired one to four spikes and then ceased firing. In Fig. 2B, a computational model demonstrating comparable firing patterns to experimental measurements. The integrate-and-fire model representation incorporated a range of M-current conductances (g M ) in discrete cells (see Table 1), which allowed for reproduction of the full range of experimentally observed action potential firing dynamics.
The distribution of neuronal firing resulting from simulation in the ANS layered network model is shown in Fig. 3. In Fig. 3A, the blue rectangles represent the sympathetic branch (SNS) of the autonomic nervous system with central nervous system input (CNS), the intrathoracic network layer (ITNS), and the intracardiac nervous system layer (ICNS). The parasympathetic branch (PNS) is shown in green with a CNS input layer connected to the ICNS layer, representing the direct vagal connection. Each layer contains a random network of model generated integrate-and-fire neurons with delayed rectifier potassium, leakage, M-type potassium and synaptic currents. Synaptic currents constitute and generate the intra-and inter-network connections. Inter-network connections are represented by arrows. Within each box, the predicted distribution of firing types from the model network is shown and was noted after the system reaches steady-state firing, averaged over seven 60-s simulations. The distribution of neuronal firing types is classified empirically by comparing firing probability over time relative to the cardiac phase relations (onset of left ventricular pressure) as shown in Fig. 3B. Model-generated neuronal firing patterns that did not clearly fit into one classification are labelled as 'other' . Figure 3B shows examples of neural firing relative to the phase of the cardiac cycle in ICNS neurons recorded from dogs (Beaumont et al., 2013). Temporal correlation to left ventricular pressure was used to classify firing types as follower, phasic or tonic.
We next set out to model and simulate the impact of sympathetic and parasympathetic stimulation on heart rate shown in Fig. 4 left panels. We employed the Behar-Yaniv model of the rabbit SAN as the model incorporates both adrenergic and muscarinic receptor pathways and subcellular signalling events associated with activation of the pathways (Behar et al., 2016). In Fig. 4A, the simulated effect of application of SNS stimulation (time course shown in the middle panel in blue) resulted in an expected increase in heart rate (pink trace in middle panel) over 60 s of stimulation. The predicted increase in heart rate shows good agreement when compared to experimental data (black symbols) from  (n = 4).
The bar graph in the right panel shows that simulated peak heart rate during SNS stimulation agrees with experimental data from (Ng et al., 2001). In Fig. 4B middle panel, the application of PNS stimulation (green) resulted in prediction of a rapid reduction in heart rate (pink trace) during 60 s stimulation. The simulated minimum heart rate during PNS stimulation was comparable to experimental data from (Ng et al., 2001) (right panel). In Fig. 4C, the application of SNS stimulation followed by PNS stimulation for 60 s in the computational model resulted in a rapid reduction in heart rate upon PNS stimulation, indicating the dominance of the parasympathetic branch to control heart rate. Notably, cessation of PNS stimulation results in immediate heart rate increase from continued stimulation of the sympathetic branch of the ANS. The firing dynamics of SAN action potentials following removal of PNS stimulation are shown in the right panel.
In Fig. 5 we show the computational model representation and predicted effects of sympathetic and parasympathetic nerve stimulation (SNS and PNS) on ventricular action potential duration (APD 80 ) and Ca 2+ transient (CaT) in the Soltis-Saucerman rabbit model under conditions of constant pacing. To study J Physiol 601.17 the impact of the ANS inputs only on the ventricular myocytes (without a change in the pacing frequency), we set pacing cycle length to 320 ms for all simulations. In Fig. 5A, the effect of SNS stimulation at a pacing cycle length of 320 ms shows a predicted reduction in the APD 80 due to K channel modification (Chen et al., 2014; and marked increase in CaT over the duration of the stimulation with increased maximum and minimum calcium transient amplitude (green lines in right panels) during the whole stimulation range. The CaT increase is due to β-adrenergic stimulation enhanced intracellular Ca 2+ uptake and release. Our simulations show a reduction in the APD 80 and CaT when the SNS stimulation ends and then a return to close to baseline. In Fig. 5B, the contrasting effect of PNS stimulation at a pacing cycle length of 320 ms is shown. The model predicts that both APD 80 (middle panel) and maximum and minimum calcium transient amplitude (green lines in right panels) with PNS stimulation alone was unchanged compared to baseline. Finally, in Fig. 5C, we show the predicted effect of constant SNS stimulation, and then application and withdrawal of PNS stimulation after 20 s (during concomitant SNS application). The predicted APD 80 and CaT were more complex in this situation, with an initial decrease in APD 80 following application of SNS stimulation, and then an increase due to transient application of PNS stimulation.
Following the development of the models of sympathetic and parasympathetic branches and the effects of stimulating these branches on both the rabbit SAN alone and the rabbit ventricular model alone, we next set out to combine these models and allow the Behar-Yaniv SAN model to couple to and drive pacing in the Iancu-Soltis-Saucerman ventricular model that is also subject to additional cell-based effects of ANS stimulation (Fig. 6). In Fig. 6A, the first combined model representing the impact of SNS stimulation on the SAN and the resulting SAN pacing frequency used to drive ventricular pacing is shown. The time course of the effects of SNS stimulation is shown in the blue trace in the middle panels. The simulated effect of dynamic SAN pacing combined with direct effects of simulated 60 s SNS stimulation on the ventricular APD 80 is shown in the red trace in the middle panel. Experimental data from the Ripplinger Lab  is also shown in black symbols and demonstrate good agreement with the predicted outputs from the coupled model. The maximum and minimum amplitude (green lines) of calcium transients during the whole stimulation range are shown in the right panels.  4). Right panel shows that simulated peak heart rate during SNS stimulation agrees with experimental data from Ng et al. (2001). B, PNS stimulation (top panel) was applied for 60 s, and heart rate was predicted to decrease to ∼140 bpm. Simulated minimum heart rate during PNS stimulation comparable to experimental data from Ng et al. (2001) (right panel). C, simulation shows effect on heart rate when SNS stimulation was applied continuously and PNS stimulation was applied transiently for 60 s (top panel in C). Heart rate was reduced to ∼200 bpm during PNS and then recovered following removal of PNS. AP firing dynamics during removal of PNS are shown in the right panel. [Colour figure can be viewed at wileyonlinelibrary.com] In Fig. 6B, a simulated 60 s stimulation of the parasympathetic branch of the ANS and its impact on both the SAN pacing frequency and the pacing induced changes combined with direct PNS effects on ventricular APD 80 is shown in red. The simulations demonstrate a small increase in APD 80 but also a notable effect on CaT (in contrast to PNS effect on a constantly paced ventricular tissue shown in Fig. 5B) as discussed above. This suggests that slower pacing from the PNS effect on the SAN leads to a notable increase in ventricular action potential duration and an effect on calcium transient amplitude. Again, the maximum and minimum amplitude (green lines) of calcium transients during the whole stimulation range are shown in the right panels. In contrast to SNS-induced increase in the Ca 2+ transient, PNS stimulation causes an abrupt decrease in Ca 2+ transient amplitude due to I CaL reduction. Our simulations also show that PNS stimulation results in reduction in I Ks and I Kr due to inhibition of phosphorylation and slower pacing frequency (reverse rate dependence) as shown in Fig. 6D.
Finally, in Fig. 6C, SNS stimulation is applied to the SAN resulting in a change in heart rate that drives the pacing frequency in the ventricular cells. During SNS stimulation, a transient PNS stimulus (60 s) is applied as shown in green. The added effect of PNS activation is to slow SAN pacing frequency and consequently the frequency of pacing of the ventricular cells. The pacing effect is combined with and interacts with subcellular signalling events that occur as a result of postsynaptic receptor activation in the ventricular cells. Shown in the middle are predictions where SNS stimulation was applied through the whole simulation, and PNS stimulation was transiently applied between 20 s and 80 s. Model predictions indicate that APD 80 was decreased following SNS stimulation, and then partially recovered with PNS stimulation. Maximum and minimum amplitude (green lines) of calcium transients during the whole stimulation range are shown in the right panel. During the SNS stimulation phase that induces a faster pacing rate, a shorter APD ensues as well as marked functional upregulation of the L-type Ca 2+ channel from phosphorylation, which leads to enhance Ca 2+ transient amplitude and increase SR Ca 2+ -ATPase (SERCA) activity. On the other hand, PNS stimulation was shown to slow rate, increase APD and block the phosphorylation-mediated changes in intracellular calcium. However, the interplay between the complex signalling events and rate dependence will need further

Figure 5. Simulated effects of sympathetic and parasympathetic stimulation on action potential duration (APD 80 ) and Ca 2+ transient (CaT) in the rabbit ventricular computational myocyte model during periodic constant pacing
Maximum and minimum amplitude (green lines in right panels) of calcium transients during the whole stimulation range are shown for each case. A, SNS stimulation was applied at cycle length of 320 ms, and APD 80 was decreased. B, PNS stimulation was applied at cycle length of 320 ms. APD 80 was unchanged from baseline. C, simulated SNS stimulation was applied throughout the simulation, while PNS stimulation was transiently applied after 20 s. The model was paced at cycle length of 320 ms. APD 80 first decreased due to application of SNS stimulation and then increased after addition of PNS stimulation.
[Colour figure can be viewed at wileyonlinelibrary.com] J Physiol 601.17 exploration in future studies to identify the parameter spaces that define 'imbalance' that may be favourable for arrhythmia formation.
In Fig. 7, we demonstrate a multiscale model utilizing the 'atom to the rhythm' approach in our recent study (Yang et al., 2020). In this case, molecular dynamics simulations were undertaken using an atomistic structural model of the β-adrenergic receptor (βAR) in order to predict quantitative estimates of noradrenaline (NA) binding affinity and association/dissociation ('on'/ 'off ') rates to the βAR (Fig. 7A, left). These molecular dynamics simulation-derived rates provided us with a reasonable estimate of the time course of β-adrenergic receptor-noradrenaline binding, with ∼97% saturation in a 100 s range for a nearly saturating ligand concentration of 10 μM. Interestingly, using recently published experimental 'on' and 'off ' rates for noradrenaline-β-adrenergic receptor interactions (Xu et al., 2021) resulted in substantially slower kinetics and only ∼17% receptor saturation at 10 μM ligand concentration in ∼1 s and nearly complete receptor saturation only on an ∼1 h time scale. We did not see appreciable receptor binding for lower noradrenaline concentrations on a second-long simulation time using that model. On the contrary, using our molecular dynamics simulation-derived model we observed gradual decrease of receptor saturation with a longer time course for lower concentrations as expected.
The molecular dynamics simulation derived rates were then used to populate the function scale model of NA-βAR interaction in the combined Iancu-Soltis-Saucerman ventricular cardiomyocyte model (see Table 2). This new function scale model was then used to make predictions of electrical activity in a one-dimensional cardiac ventricular tissue (1.65 cm in length) in the absence of autonomic input (baseline) as shown in the computed pseudo-electrocardiograms (pseudo-ECGs) from the tissue in Fig. 7B (middle). Pacing frequency was set to heart rate generated by the Behar-Yaniv SAN model. We then applied a different protocol in Fig. 7B (right), showing the effect of SNS stimulation (blue) on the computed pseudo-ECG followed by addition of PNS stimulation (green), which was applied in the model simulation between 20 and 80 s (SNS was applied for the whole simulated time course of 100 s). Again, the changes in pacing frequency are driven by the autonomically mediated changes in the SAN model. The red line indicates the increased amplitude of the T wave during PNS stimulation. Computed electrograms at different time points are shown at the bottom with the peak of the T wave indicated by the red circles. These simulations indicated that T wave amplitude increases during PNS and returns after withdrawal.
In Fig. 8, we predict the effect of cessation of SAN input to the ventricular model in the merged multiscale Iancu-Soltis-Saucerman model. Pacing frequency was set to the heart rate generated by SAN model until 415 s and then stopped to simulate a sinus pause. In the absence of SNS application and cessation of beating at 415 s, the simulated non-diseased (Fig. 8A) and diseased hearts (Fig. 8B) are predicted to have expected normal behaviour with no additional electrical or Ca 2+ activity following decoupling of the SAN from the ventricular model representation (top row). However, in response to a sympathetic surge as shown in Fig. 8C, even a non-diseased heart generates a single spontaneous AP (red peak in middle) that is triggered by a delayed afterdepolarization, following SNS surge and cessation of 10% decrease Human Bossuyt et al. (2005), Bundgaard & Kjeldsen (1996), Shamraj et al. (1993) beating at 415 s. In Fig. 8D simulations in a diseased heart (simulated heart failure) with the same protocol as in Fig. 8A predict the emergence of multiple triggered afterdepolarizations. In Fig. 9, we predict the effects of adding parasympathetic nervous system stimulation to the proarrhythmic effects of sympathetic surge shown under the same conditions as in Fig. 8. In the setting of overly active sympathetic stimulation in non-diseased heart (blue bar at top of Fig. 9A) via the activation of βAR, we applied PNS stimulation (green bar in top of Fig. 9A). The dynamics of the single-cell AP and CaT is shown. In Fig. 9A PNS stimulation was applied and eliminates the triggered beat in a non-diseased heart that was predicted in Fig. 8. In Fig. 9B a simulation in a diseased heart (simulated heart failure) with the same protocol as in Fig. 9A resulted in suppression of triggered activity at 422 s compared to continued firing in Fig. 8 with no PNS stimulation.
Finally, in Fig. 10, we employed the atom to tissue model and tested the case in 1D tissue that included the predicted rates of NA interaction with βAR from our atomistic structural model and simulated the effect of the PNS stimulation on spontaneous beat triggering events. The simulations suggest the triggered activity was initiated and maintained for a much longer period of time in the diseased heart -for about 20 s (Fig. 10B) compared to the observation of a single spontaneous beat in the healthy heart (Fig. 10A). When PNS stimulation was applied back in the model after termination of pacing, the triggered activity in the diseased heart was discontinued within 5 s as shown in Fig. 10D and was not initiated in the healthy heart condition (Fig. 10C).

Discussion
The ultimate purpose of our study is to provide a digital twin for modelling, simulation and prediction of neuromodulation to treat, reverse and prevent heart disease. As autonomic modulation is increasingly shown to be a powerful tool for cardiac arrhythmia prevention and therapy, new tools to allow finer grained prediction of appropriate interventions are needed and the tool we present represents a critical step in computational modelling and simulation (Dacey et al., 2022;Dusi et al., 2021;Hadaya et al., 2022;Hanna et al., 2018;Kember et al., 2017;Liu et al., 2018;Zhu et al., 2022).
In this study, we set out to advance prediction of the effects of stimulation of the autonomic nervous system on cardiac electrophysiology. To do so, we developed a multiscale digital twin model to capture subcellular, cellular and system scale features of autonomic control of the heart. Some key features of the new modelling and simulation approach include a multi-layered structure of the autonomic nervous system that represents both the sympathetic and the parasympathetic branches as shown in Figs. 1 and 3. The layers include sparse random intralayer connectivity, synaptic dynamics and conductance-based integrate-and-fire cellular dynamics in close agreement with experiment (Springer et al., 2015) as shown in Fig. 2.
Another key feature of the neurocardiac computational model is the modelled connection between the autonomic nervous system and both the cardiac pacemaker and the contractile cells as shown in Figs. 4-6. The result is a dynamic model where modification of pacemaker frequency drives initiation of electrical signals in the contractile cells. Receptor-mediated signalling on both pacemaker and contractile cells leads to a cascade of subcellular cardiac signalling processes occurring on multiple time scales that result in predicted functional changes to cardiac electrophysiology in both nodal and contractile cell representations.
We show the predicted effect of application of SNS stimulation resulted in an expected increase in heart rate (pink trace in middle panel) over 60 s of stimulation in Fig. 4. The predicted increase in heart rate shows good agreement when compared to experimental data (black symbols) from  (n = 4), but the experimental data show a decrease in the action potential duration following the peak. There are several possible mechanisms of the heart rate reduction following the peak increase in the experiment including β 1 AR desensitization, depletion of nerve terminals due to persistent high frequency stimulation, and/or washout of neurotransmitters and precursors during the experiment. We have recently developed and published a reduced yet detailed computational model of the β 1 -adrenergic signalling cascade to a system of two differential equations by eliminating extraneous variables and applying quasi-steady state approximation . The structure of the reduced model reveals that the large cAMP transients associated with abrupt β 1 AR activation are generated by the interplay of production/degradation of cAMP and desensitization/resensitization of β 1 ARs. The explicit modelling of desensitization of β 1 ARs by PKA is a likely contributing mechanism to the reduction in rate observed in the experimental data. This model is available as a component that can be used in the simulator, but we have not included it in this study.
We demonstrated that we can construct a multimodal model that includes a layered autonomic neural representation connected to nodal and contractile cells with cell signalling dynamics (Figs. 4-6). However, a major challenge of extending the model to other receptor-mediated signalling events is experimental determination of needed parameters for the model. To begin to overcome this challenge, we set out to determine if we could utilize atomic-scale approaches to generate some known model parameters that have already been constrained experimentally to allow model validation. A successful example is shown in Fig. 7, where we present an extension of the multiscale model from the atom to the rhythm. Such models were previously used in our studies to predict emergent pro-arrhythmia proclivities of small-molecule drugs, blocking cardiac potassium channel encoded by human ether-à-go-go-related gene (hERG) (DeMarco et al., 2021;Yang et al., 2020). There we used all-atom molecular dynamics simulations to predict Figure 10. Effects of transient exposure to PNS and SNS surge in simulated healthy (non-diseased) and diseased (simulated heart failure) one-dimensional cardiac tissues using a multiscale model that includes autonomic nerve stimulation on rabbit sinoatrial node (SAN) and ventricular cell models Pacing frequency was set to heart rate generated by the Behar-Yaniv SAN model, and we used predicted rates of βAR-NA interactions obtained from atomic simulations in the functional scale models. One-dimensional tissue (1.65 cm) with SNS was simulated for 1529 beats. For each set of panels, pseudo-ECGs are at the top, voltage time course and calcium transient between t = 415 s and 420 s are at the middle and bottom in blue. A, in the setting of overly active SNS in a healthy heart model, we ceased application of pacing stimuli at 415 s (the 1529th beat). We observed a spontaneous beat (red) that was triggered after cessation of pacing at t = 415 s. B, following the same protocol as in panel A, for diseased heart model, the extensive triggered activity was initiated. C, with addition of PNS stimulation applied in the non-diseased model after t = 415 s, the triggered beat was suppressed. D, with addition of PNS stimulation, the triggered activity (red) terminated around 418 s in the diseased model. Grey arrows indicate triggered action potentials and Ca 2+ transients on space-time representations of one-dimensional tissues. [Colour figure can be viewed at wileyonlinelibrary.com] association and dissociation ('on' and 'off ') rates of hERG blocking drugs to the open channel pore. They were computed based on free energy and diffusion coefficient profiles from umbrella sampling molecular dynamics simulations. Those rates were used as parameters of functional kinetic models of state-dependent hERG channel-drug interactions, which were used to reproduce experimental dose-response curves and incorporated into cardiac cell and tissue models to predict emergent pro-arrhythmia outcomes (DeMarco et al., 2021;Yang et al., 2020).
In the model pipeline shown in Fig. 7, we utilized atomic-scale molecular dynamics simulations to predict the association and dissociation ('on' and 'off ') rates of noradrenaline interactions with the β-adrenergic receptor. Since ligand binding to the adrenergic receptor and other G protein-coupled receptors does not follow a linear path (Provasi et al., 2009;Schneider et al., 2015) as a drug binding to the hERG channel pore (DeMarco et al., 2021;Yang et al., 2020), we used a different enhanced sampling molecular dynamics simulation approach, well-tempered metadynamics, and estimated dissociation rate using Kramer rate theory formalism based on computed free energy profile and diffusion coefficient at the free energy barrier computed from a separate restrained molecular dynamics simulation. Calculation of diffusion coefficient profiles across an entire reaction coordinate as was done in previous studies (Berezhkovskii et al., 2011;DeMarco et al., 2018;DeMarco et al., 2021;Setny et al., 2013;Vorobyov et al., 2014;Yang et al., 2020;Zhu & Hummer, 2010;Zhu & Hummer, 2012) can help refine computed ligand-receptor 'on' and 'off ' rates and will be explored in our subsequent studies along with other rate computation methods (Meral et al., 2018;Palacio-Rodriguez et al., 2022;Pang & Zhou, 2017;Wang et al., 2023). This methodological shortcoming may contribute to a substantial difference between our MD estimated and substantially slower experimental 'on' and 'off ' rates from a recent study (Xu et al., 2021). Achieving agreement between computed and experimental protein-ligand rates has proven to be a challenging task with results varying by several orders of magnitudes and discrepancies attributed to issues with enhanced sampling methodologies, choice of reaction coordinate (s), force field accuracy, presence of hidden barriers and other MD simulation related factors (Wang et al., 2023) as well as those influencing accuracy of experimental rate measurements using, e.g. competitive binding assay (Georgi et al., 2019;Hoare, 2021;Sykes et al., 2019).
Despite these challenges, the molecular dynamics simulation-derived rates predicted a nearly identical time course of cAMP production compared to the original Soltis-Saucerman and Iancu-Harvey subcellular signalling cardiac myocyte models, where in the latter quasi-steady-state approximation (instantaneous binding of NA to the receptor) was assumed. This indicates that βAR-NA interaction is not a rate-limiting step for cAMP production in our molecular dynamics informed model. This is also in agreement with our recent model sensitivity study , which indicated that production/degradation of cAMP and desensitization/resensitization of βARs are rate-limiting steps for adrenergic stimulation.
Thus, utilizing these predicted rates resulted in what appears to be an effective approach to constrain model parameters through simulation. Indeed, our predicted parameters yielded nearly identical outputs in the multi-scale model compared to those constrained by experimental data alone as shown in Figs. 4-6. This gives us some confidence that the approach can be extended to consider multiple conformational states and other types of G protein-coupled receptors in future studies and to potentially pursue quantitative investigations into adrenergic receptor binding of different ligands. For example, future studies might allow for comparison between biased agonism versus biased antagonism using a molecular dynamics informed approach similar to the one we use here.
In healthy individuals, the autonomic nervous system plays a critical role in regulating cardiac function, with the PNS primarily responsible for controlling heart rate and the SNS regulating both heart rate and contractility. However, chronic cardiac injury can lead to an exaggerated sympathetic response. In Figs. 8-10, we explored the potential for activation of the parasympathetic nervous system to suppress arrhythmia triggers induced by an overactive sympathetic nervous system. The simulation is derived from earlier studies suggesting that high sympathetic tone can be combined with muscarinic activation to reduce triggered activity in cardiac tissue (Harvey & Belevych, 2003;Mehra et al., 2022;Meijborg et al., 2020;Silvani et al., 2016). The emergence of arrhythmogenic triggered beats observed in the presence of maximal SNS stimulation was more prominent in a diseased heart model mimicking heart failure condition as shown in Fig. 8B. Although we are not able to confirm maximal physiological range of diastolic Ca 2+ , it is generally accepted that maximal adrenergic stimulation is likely to be a key mechanism contributing to Ca 2+ overload conditions as our study suggests. Notably, we also observed a protective role of reapplied parasympathetic stimulation, which substantially reduces triggered beat activity in the diseased heart and eliminates it altogether in the healthy heart as shown in Figs. 9 and 10.
Such simulations may be useful to extend and improve understanding underlying molecular mechanisms of vagal nerve stimulation therapies, which are emerging techniques for non-invasive and non-pharmacological arrhythmia treatments (De Ferrari & Schwartz, 2011;J Physiol 601.17 Hanna et al., 2018;Liu et al., 2020;Zhang & Mazgalev, 2011). The model parameters can be tuned to be applied in real cases in the absence of drug therapy, or more likely, in combination with therapeutic pharmacology. However, it is important to note that chronic cardiac injury can lead to alterations in autonomic regulation of the heart, including a loss of central vagal drive, which may also require consideration (Francis Stuart et al., 2018). Overall, incorporating the effects of chronic cardiac injury on the autonomic nervous system can be challenging and requires careful consideration of the underlying mechanisms and their potential interactions with experimental conditions.
The present results suggest that antagonism of SNS responses occurring at the myocyte level is one factor that may contribute to the potential for PNS stimulation to attenuate arrhythmogenic responses. However, the interactions between the PNS and SNS are more complex than the steady-state responses reflected in the simulations presented here. For example, dynamic changes in SNS and PNS tone that occur in disease states such as sleep apnoea and epilepsy can produce complex temporal responses that may contribute to the generation of triggered activity (Chadda et al., 2018;Costagliola et al., 2021;Geovanini & Lorenzi-Filho, 2018;May et al., 2017;Soh et al., 2022). Furthermore, the current models do not account for the changes in presynaptic or postsynaptic signalling mechanisms that are known to occur with myocardial infarction and heart failure. One fundamental experimental and simulation study from Iancu and Harvey demonstrated complex functional responses to activation of the autonomic nervous system arising from the interplay of signalling cascades resulting from stimulation of β 1 -adrenergic receptors (β 1 AR) and M2 muscarinic acetylcholine receptors (M 2 R) and resulting compartmentalized cAMP signalling in adult cardiac myocytes (Iancu et al., 2007(Iancu et al., , 2008. While it is well known that the M 2 R plays an important role in parasympathetic regulation of cardiac myocyte function by modulating cAMP production (Harvey & Belevych, 2003) through G i -dependent inhibition of adenylyl cyclase (AC) activity, Harvey's group showed the effect muscarinic stimulation has on ventricular myocytes is more complex because M 2 R stimulation can also activate a stimulatory pathway through stimulation of AC4/7 (Sunahara et al., 1996;Tang & Gilman, 1992;Taussig & Gilman, 1995). Modelling and simulation showed that representing cAMP production in different microdomains with different kinetics can account for the experimental observations. We have included these domains in the model presented here to show how complex subcellular signalling manifests during autonomic imbalance.
Furthermore, we still need to learn how to fine-tune vagal nerve stimulation to provide safe and efficacious patient-specific treatment of cardiac arrhythmias resulting from the autonomic imbalance associated with various disease states. This would require an accurate prediction of time-dependent and autonomic nerve stimulus-specific responses on heart rate and stroke volume, for which our developed neurocardiac simulator may provide a first framework that can be easily extended and improved to include a variety of new details relevant to new questions (Mehra et al., 2022;Sridharan et al., 2022).
A very important element that will be included in the future state is the remodelling of the ANS that occurs during a variety of disease states. While we did not address all of these issues in the coarse-grained approach we took to 'disease' in this first study, we do recognize its importance as a current limitation and also an important future direction. One important point to emphasize is that the digital twin can be utilized by anyone in the public domain and so critical questions such as how does disease-induced remodelling of the nerves modify the predicted effect of neuromodulation? Incorporating new experimental work by Habecker and others on nerve remodelling will allow the neurocardiac digital twin to predict expected responses to autonomic inputs during disease states based on relevant data (Clyburn, Andresen et al., 2022;Clyburn, Sepe et al., 2022;Mehra et al., 2022;Sepe et al., 2022). It is important to note that the exact mechanisms underlying some disease-induced nerve remodelling and other changes are not fully understood but can be continuously incorporated as new data become available. A full range of parameter regimes can be explored through the digital twin, including the potential effects of changing discharge frequencies, neurotransmitter concentrations and neuropeptide release patterns that are altered during various disease progressions.
In terms of models of cardiac function, there are several approaches that can be taken to incorporate these changes. One approach is to use experimental data to parameterize models that account for altered autonomic input. For example, models can be adjusted to simulate the effects of reduced PNS input on heart rate and contractility, or increased SNS input on contractility. Another approach is to use computational models that incorporate physiological data to simulate the effects of chronic cardiac injury on autonomic regulation. For example, these models can take into account changes in baroreflex sensitivity, which can lead to an exaggerated sympathetic response to changes in blood pressure. It is important to note that the specific approach taken to incorporate changes in autonomic regulation into models of cardiac function will depend on the research question being addressed and the available data. Overall, however, understanding the effects of chronic cardiac injury on autonomic regulation of the heart is an important area of research, as it has important implications for the development of new therapies for cardiovascular disease. Indeed, these studies comprise the future applications of the digital twin for neurocardiac modulation, which can be modified to include specific elements of an individual or group disease state.
In summary, the developed neurocardiac simulator allows us to directly predict the response of ANS stimulation on heart rhythm through coupled SAN and ventricular cardiac tissue effects in both healthy and disease conditions such as heart failure. It can be used to predict, prevent and reverse resultant cardiac arrhythmia through precisely targeted parasympathetic and/or sympathetic nervous system stimulations, which is a major goal of our collaborative mutli-scale project supported by NIH Common Fund SPARC initiative. In the future, translation of experimental data on ANS stimulation in animal models to human cardiac electrophysiology (Morotti et al., 2021) can be applied to supplement scarcer human experimental data. Moreover, the digital twin can be personalized through use of patient-specific induced pluriopotent stem cell derived cardiomyocytes (iPSC-CMs) and translation of their immature phenotype to mature cardiac myocyte electrophysiology and calcium dynamics as we showed recently (Aghasafari et al., 2021).