Modeling effects of voltage dependent properties of the cardiac muscarinic receptor on human sinus node function

The cardiac muscarinic receptor (M2R) regulates heart rate, in part, by modulating the acetylcholine (ACh) activated K+ current IK,ACh through dissociation of G-proteins, that in turn activate KACh channels. Recently, M2Rs were noted to exhibit intrinsic voltage sensitivity, i.e. their affinity for ligands varies in a voltage dependent manner. The voltage sensitivity of M2R implies that the affinity for ACh (and thus the ACh effect) varies throughout the time course of a cardiac electrical cycle. The aim of this study was to investigate the contribution of M2R voltage sensitivity to the rate and shape of the human sinus node action potentials in physiological and pathophysiological conditions. We developed a Markovian model of the IK,ACh modulation by voltage and integrated it into a computational model of human sinus node. We performed simulations with the integrated model varying ACh concentration and voltage sensitivity. Low ACh exerted a larger effect on IK,ACh at hyperpolarized versus depolarized membrane voltages. This led to a slowing of the pacemaker rate due to an attenuated slope of phase 4 depolarization with only marginal effect on action potential duration and amplitude. We also simulated the theoretical effects of genetic variants that alter the voltage sensitivity of M2R. Modest negative shifts in voltage sensitivity, predicted to increase the affinity of the receptor for ACh, slowed the rate of phase 4 depolarization and slowed heart rate, while modest positive shifts increased heart rate. These simulations support our hypothesis that altered M2R voltage sensitivity contributes to disease and provide a novel mechanistic foundation to study clinical disorders such as atrial fibrillation and inappropriate sinus tachycardia.


Introduction
The cardiac muscarinic receptor (M2R) plays a crucial role in regulating heart rate variability and vulnerability to atrial arrhythmia by modulating the acetylcholine (ACh) activated K + current I K,ACh . Cardiac K ACh channels are heteromultimers composed of two G-protein-coupled inward rectifier K + channel subunits, Kir 3.1 and Kir 3.4 [1]. ACh activation of M2R triggers dissociation of the G beta-gamma subunits (G βɣ ) that in turn directly activate Kir 3.x subunits to conduct I K,ACh . Unexpectedly, M2Rs were discovered to possess an intrinsic ability to sense transmembrane voltage [2] and the affinity of the receptor for ligands was noted to vary in response to changes in membrane voltage [3]. In particular, the affinity of the receptor for ACh is increased at hyperpolarized membrane potentials and decreased at depolarized potentials. The changes in affinity exert a downstream effect on the K ACh channel such that the channel is more active (more current) at hyperpolarized potentials and less active (less current) at depolarized potentials. The observation that M2Rs are intrinsically voltage sensitive has profound implications for cellular signaling in excitable tissues, such as heart. For example, voltage sensitive behavior provides a mechanistic explanation for a decades-old enigmatic process called I K,ACh "relaxation" gating. Relaxation gating refers to a time-dependent change in current magnitude following a depolarizing or hyperpolarizing voltage step [4] and has important consequences for shaping the cardiac action potentials (AP), especially in the sinus node. We recently proposed that relaxation gating represents a voltage dependent change in ACh affinity induced by voltage dependent conformational changes within M2R [5]. Our experimental data provide a mechanistic basis to explain the participation of I K,ACh in the modest chronotropic effects induced by resting vagal tone. As a result of conformational changes in the M2R, the affinity for ACh varies throughout the cardiac electrical cycle such that low (subsaturating) ACh concentrations preferentially activate I K,ACh during diastolic membrane voltages thereby slowing the spontaneous firing rate without appreciably altering AP duration (APD).
Alterations in the voltage sensitivity of M2R could theoretically contribute to cardiovascular diseases that clinically present with apparent changes in vagal tone. For example, genetic variants in M2R that shift the receptor occupancy into the hyperpolarized state would be expected to increase the affinity of the receptor for ACh and thus activate more K ACh channels at a given ACh concentration (or degree of vagal tone). Accordingly, genetic variants in M2R that shift the receptor occupancy into the hyperpolarized state might explain the clinical phenotype of vagally-mediated atrial fibrillation (AF), patients who present with bradycardia in the setting of physiological (basal) ACh concentrations. Alternatively, genetic variants in M2R that shift the receptor occupancy into the depolarized state would be expected to decrease the affinity of the receptor for ACh and thus fewer K ACh channels activate at a given ACh concentration (or given degree of vagal tone). This would decrease the effects of vagal modulation of heart rate, thereby increasing basal heart rate, as observed in the syndrome of inappropriate sinus tachycardia (IST).
To provide insights into the contribution of M2R voltage sensitivity to cardiac electrophysiology in physiological and pathophysiological conditions, we extended our previous Markovian model of M2R [5] to incorporate G βɣ -mediated activation of the K ACh channel and integrated the revised Markovian model into a human model of the sinus node (SN) cell model [6]. Based on experimental data from isolated human SN cells [7,8], Fabbri and colleagues recently published a comprehensive model of the human SN pacemaker cell that faithfully recapitulated the effects of autonomic modulation as well as mutations associated with SN dysfunction [6]. In the Fabbri model and its parent model [9], I K,ACh is described by a voltage-and [ACh]-dependent gate, but the intrinsic voltage sensitivity of M2R is not incorporated. Here, we introduce a model reproducing the effects of M2R voltage sensitivity on human SN cell function under physiological and pathophysiological conditions. These simulations support our hypothesis that altered M2R voltage sensitivity contributes to disease and provide a novel mechanistic foundation to study clinical disorders such as AF and IST.

Effects of low concentrations of ACh on sinus node action potentials
For decades, the contribution of I K,ACh to the modest chronotropic effects of 'physiological' or low-dose ACh has been debated [10][11][12]. Based on the M2R voltage-dependent properties, we predict that subsaturating ACh concentrations exert a larger effect during diastolic (hyperpolarized) membrane voltages, compared to the voltages during the cardiac AP (depolarized) [3,5]; thus preferentially slowing the pacemaker rate with minimal effect on APD. To test this hypothesis, we simulated the effects of varying low, sub-saturating concentrations (e.g., 20-100 nM) of ACh on sinus node AP properties (Fig 1). The most striking effect of increasing ACh concentration was reduced slope of phase 4 depolarization and the corresponding increase in basic cycle length ( Fig 1A, Table 1, V shift = 0), with minimal shortening of APD 90 ( Table 2, V shift = 0). Thus, the basic cycle length (BCL) increased from 827 ms in the absence of ACh, to 1585 ms in the presence of 0.1 μM ACh. The AP amplitude decreased steadily by up to 4.7 mV at the highest tested concentration of 0.1 μM ACh (inset of Fig 1A). Simulated open probability O and I K,ACh for different ACh concentrations are shown in Fig 1B and 1C. Taken together, these simulations indicate that subsaturating concentrations of ACh slow spontaneous excitation of the SN cell by inhibiting the rise of phase 4 depolarization, without appreciably shortening APD 90 or reducing the amplitude of the AP.

Modeling perturbations in M2R voltage sensitivity
In our previous experimental studies using isolated feline left atrial myocytes, the voltage dependence of M2R was explored by measuring the ACh concentration-I K,ACh response relationship at hyperpolarized (-100 mV) and depolarized membrane voltages (+50 mV) [3,5]. These experiments indicated that the affinity of the receptor for ACh was greater at hyperpolarized membrane voltages, compared to depolarized voltages. We reasoned that, similar to voltage-gated ion channels, putative disease-associated mutations in M2R might alter the voltage sensitivity of the M2R, with unique consequences for sinus node AP properties and heart rate responses. We modified rate parameters (Eqs 5 and 6) to shift the receptor occupancy towards a hyperpolarized state (higher affinity) or a depolarized state (lower affinity) (S1A Fig). Thus, we simulated the effects of shifting the M2R voltage sensitivity (Fig 2, S3C Fig). S1B Fig highlights the effects of voltage shifts on the state O. Negative shifts (e.g., V shift of -30 mV and -150 mV), which caused a more hyperpolarized state of the receptor, increased the occupancy of the U1 and B1 states, as well as the state O. Negative voltage shifts resulted in a slight leftward shift in the concentration-response curve when the cell was held at +50 mV ( Fig 3E). Likewise, positive shifts in M2R voltage sensitive parameters caused a slight rightward shift in the concentration-response curve for a holding potential of V h = -100 mV (Fig 3A). To avoid local minima during the parameter fitting, the model was forced to favor the open state (state O equal to 1), at the maximum concentration of 10 μM ACh and to favor the closed state (O equal to 0), with no ACh present. Further, U1 was forced to be as high as possible at a holding potential of -100mV and U2 at +50 mV, in the absence of ACh. Thus, negative and positive voltage shifts did not move the steady-state concentration-response relationships outside of these ranges. Notwithstanding, the voltages experienced by the single cell model vary between -60 and +30 mV (Fig 1A), within the minimum and maximum ranges defined by the parameter fitting. We previously described the kinetics of "relaxation" gating of I K,ACh in terms of activation and deactivation of K ACh channels in the setting of subsaturating ACh concentrations [5]. Activation kinetics were measured by first stepping to a depolarized voltage (+60 mV) to close a large portion of K ACh channels at a physiological voltage, followed by stepping through a range of voltages to measure activation of I K,ACh . Deactivation kinetics were assessed by a prepulse to a hyperpolarized voltage (-100 mV) to open channels, followed by variable test voltage steps to measure the rate of K ACh channel closure. Accordingly, simulated I K,ACh evoked by 0.1 μM ACh using the activation and deactivation voltage protocols are presented in S1 Fig and the kinetic parameters are described in Table 3. The simulations recapitulate the experimental features of I K,ACh relaxation gating [5], as shown in S1 Fig Next, we simulated the effects of negative and positive voltage shifts in M2R voltage-sensitive parameters on sinus node APs, together with the corresponding effects on state O and I K, ACh (Fig 4). Negative shifts (e.g., V shift of -10 mV and -30 mV), which would be predicted to increase ACh affinity, reduced the slope of phase 4 depolarization in a concentration-dependent manner, with minimal effects on AP amplitude or APD 90 (Fig 4A versus 4D, Tables 1 and  2). The reduction in the slope of phase 4 depolarization induced by a negative V shift was due to an increase in the open probability of K ACh channels relative to the control condition, thereby increasing the magnitude of I K,ACh . These results indicate that shifts in the M2R voltage sensitivity impact the rate of spontaneous depolarization of sinus node APs. To further characterize the physiological consequences of positive and negative voltage shifts, we studied the effects of variable M2R voltage sensitivity on the ACh concentrationheart rate response relationship (Fig 5A). Our model recapitulates experimental data [12] indicating that ACh concentrations ranging from 0.01 to 0.1 μM induce slowing of the spontaneous activity by 10% and 45%, respectively. Hyperpolarizing shifts in the M2R voltage dependent parameters shifted the relationship toward progressive spontaneous activity slowing. By contrast, depolarizing shifts in these parameters antagonized the spontaneous activity slowing induced by ACh. Taken together, these results suggest that shifts in M2R voltage sensitive parameters exert significant physiological effects on SN firing rate and AP parameters.  Next, we quantified the effects of the individual ACh-influenced cardiac currents on heart rate slowing. Fig 5B illustrates the relative contributions of ACh-sensitive currents, including our new model of I K,ACh , together with formulations of the L-type Ca2+ current (I Ca,L ) and the hyperpolarized-activated 'funny current' (I f ) from Fabbri et al. [6]. Inspection of the individual   ACh. (B), The individual contributions of ACh-sensitive currents to the measured concentration-heart rate response. Our current model that incorporates the voltage-sensitivity of the M2R into I K,ACh (new), together with formulations for the ACh-sensitive currents I Ca,L and I f from Fabbri et al [6] provides the best fit to the measured data. I K,ACh contributes to about 50% of the heart rate response at 0.03 μM ACh and contributes a greater proportion to heart rate slowing at higher ACh concentrations. contributions of the ACh-sensitive currents reveals that our model of I K,ACh (incorporating the voltage-sensitivity of M2R) accounts for roughly 50% of the slowing at~30nM and more for higher concentrations. These results highlight the important contribution of M2R voltage sensitivity to heart rate slowing induced by ACh.

Discussion
The observation that M2Rs are intrinsically voltage sensitive has profound implications for cellular signaling in excitable tissues, such as heart. This could have important consequences for cardiovascular drug development in that the affinity for ligands may vary in a voltage-and ligand-specific manner. The properties of M2R voltage dependence provide a novel molecular basis to explain the previously perplexing "relaxation" gating of I K,ACh , originally described in the 1970s [4]. Precisely how the receptor's voltage dependent properties influence the firing rate and shape of sinus node APs under physiological conditions, or in the presence of a M2R genetic defect, remains unknown. In order to address this question, we enhanced our previous Markovian model of M2R [5] by a 2-state Markovian model to incorporate G-protein activation of the K ACh channel. The parameters of this revised model were adapted to experimental data from isolated feline left atrial myocytes, see Methods. We integrated the revised model into a human SN single cell model and used computational simulations to gain insights into effects of voltage sensitivity of M2Rs.
While there is no controversy as to the cardiac effects of strong vagal stimulation, the participation of I K,ACh in mediating the purely chronotropic effects of low (or "physiologic") ACh concentrations has been debated over the past thirty years. Low concentrations of ACh (e.g., below 100 nM) and weak vagal stimulation reduce spontaneous pacemaker rate without altering APD or increasing the maximal diastolic potential [10]. DiFrancesco et al. proposed that weak vagal stimulation primarily inhibited I f with little to no contribution from I K,ACh . Indeed, the EC 50 for ACh modulation of I f is an order of magnitude lower than that required for I K,ACh activation [11]. However, subsequent studies have corroborated an important role for I K,ACh in mediating the chronotropic effects of weak vagal stimulation and externally applied low ACh concentrations [12]. Perhaps the strongest evidence for I K,ACh contribution to basal chronotropy comes from the Kir3.4 knock-out mouse which displays a specific deficit in heart rate variability at rest [13,14]. These studies confirm that low ACh concentrations (basal vagal tone) activate I K,ACh to slow heart rate. Our experimental and simulated data provide a mechanistic basis to explain the participation of I K,ACh in the modest chronotropic effects induced by resting vagal tone. Our experimental data predict that the affinity of M2R for ACh varies throughout the AP. The simulations performed here corroborate that subsaturating ACh concentrations preferentially open the K ACh channel during diastolic membrane voltages and thereby slow the spontaneous firing rate. As seen in Fig 5B, I K,ACh accounts for roughly half of the heart rate slowing at concentrations around 30 nM, and more for higher concentrations.
In silico modeling is valuable for understanding fundamental features of physiology and pathophysiology [15][16][17]. This is especially relevant for modeling of I K,ACh in disease states. For example, numerous disease-causing mutations in ion channel genes alter the voltage sensitivity of the channel [18,19]. In light of these observations, we simulated effects of hypothetical mutation-induced positive and negative shifts in M2R voltage sensitivity to ACh on sinus node AP properties and spontaneous firing. Modest voltage shifts (V shift of ±10 mV) exerted significant effects on the slope of phase 4 depolarization and thus on spontaneous activity responses to subsaturating ACh concentrations, such as those elicited by vagal resting tone. Moderate voltage shifts (V shift of ± 30 mV) caused more profound changes in spontaneous activity. These simulations provide a proof-of-principle for a theoretical contribution of altered M2R voltage sensitivity to cardiovascular disease states associated with changes in vagal tone. For example, parasympathetic induction and maintenance of atrial arrhythmias is a welldescribed phenomenon, first reported in 1930 [20]. Moreover, increased parasympathetic tone is an initiating factor in a subset of AF patients [21]. Our simulations suggest that genetic variants in M2R that shift the receptor occupancy towards hyperpolarized states (U1/B1) associated with increased ACh affinity of the receptor increase I KACh at a given ACh concentration (or degree of vagal tone). These genetic variants might explain, in part, the clinical phenotype of vagally-mediated AF patients who present with bradycardia in the setting of physiological (basal) ACh concentrations. Also, our simulations suggest that genetic variants in M2R that shift receptor occupancy towards depolarized states (U2/B2) associated with decreased ACh affinity of the receptor reduce I KACh at a given ACh concentration (or degree of vagal tone). Thus, positive shifts in M2R voltage sensitivity would decrease the effects of vagal modulation of heart rate, thereby increasing basal heart rate, as observed in the syndrome of IST. Indeed, a recent study confirmed decreased parasympathetic tone in IST patients [22]. While speculative, the hypothesis that altered M2R voltage sensitivity is relevant for cardiac diseases provides a novel mechanistic foundation to study disorders such as AF and IST.

Limitations of the model
There are several limitations inherent in the application of the model to describe the kinetics and behavior of I K,ACh . First, the receptor-channel model was fit and optimized to recreate the kinetics of I K,ACh occurring at a concentration of 0.1 μM ACh. While this is within the range of measured concentrations of ACh, measurements at lower concentrations could enable a more accurate reconstruction of the behavior of I K,ACh and a direct comparison to the effect of I f at such concentrations. Furthermore, we acknowledge that the description of the process of dissociation of G βɣ from the M2R to the activation of the K ACh channel is highly simplified. We used a simple channel opening description with a 2-state Markovian model neglecting different binding properties of different Kir subunits or cooperativity mechanisms in binding of G βɣ . This simplification was necessary as parameters for more complex models are not identifiable with the existent experimental data. Also, the model is not able to reproduce the distinct characteristics of the deactivation protocol time constants. The experimental data indicate that the deactivation time constant is nearly voltage-independent (Fig 3H), the model predicts an increase with membrane depolarization. Nonetheless, because the other simulated features have a very high similarity to the measured values (Fig 3B, 3C, 3F and 3G) and the time constants are in the range of less than half the length of an action potential, we believe that any error introduced does not significantly influence our findings. Additionally, dissociation or binding of G βɣ in our model does not account for changes in the process due to other influences or any binding to sites other than the K ACh channel. The SN cell model recapitulates the experimental findings that ACh inhibits I f and I Ca,L , which also contribute to slowing of spontaneous pacing rate [23]. We argue that unlike I K,ACh , the M2R voltage dependent effects do not influence I f and I Ca,L on the time scale of the AP. I f and I Ca,L inhibition are mediated by inhibition of the cAMP-dependent protein kinase A cascade that functions on a much slower time scale than the APD. By contrast, our simulations demonstrate that voltage dependent conformational changes in M2R influence I K,ACh throughout the cardiac AP, modulating both firing rate and APD. Finally, because our model does not fully integrate all the components of the autonomic nervous system, it is possible that the effects of the putative genetic mutations might be mitigated by compensatory changes in the autonomous nervous system's response to changes in heart rate.

Conclusions
The recent observation that M2Rs are intrinsically voltage sensitive has important implications for understanding the physiology and pathophysiology of parasympathetic regulation of heart rate and APD. By optimizing and integrating a new Markovian model into a human SN model, we show that low ACh concentrations preferentially slow beating rate, without shortening APD, and thereby provide additional support that I K,ACh participates in the purely chronotropic effects of basal vagal tone. Moreover, we explore the effects of altered M2R voltage sensitivity and provide a proof-of-principle foundation that altered sensitivity could result in clinical manifestations of disease states such as vagally-mediated atrial fibrillation and syndrome of inappropriate sinus tachycardia. Given the importance of parasympathetic regulation of atrial vulnerability, M2Rs represent an important therapeutic target to control or prevent atrial arrhythmias.

Model of M2R and KACh channels
We developed a Markovian model to reconstruct the behavior of K ACh channels at different ACh concentrations and varying transmembrane voltages (Fig 2). The model comprises 3 submodels: (1) A Markovian model describing the kinetics of the M2R depending on different concentrations of ACh at different voltages (M2R model), (2) a Markovian model describing the activation of the K ACh channel based on dispersion of G βɣ protein from the receptor to the channel (K ACh channel model), and (3) a model of potassium current through K ACh channels.
Markovian model of M2R. This sub-model describes biophysical properties of the M2R and comprises 4 distinct states, U1, B1, U2 and B2, whose interaction describes the affinity of M2R to changes in ACh concentration [L] and transmembrane voltage V m (Fig 2, left side). The sub-model was published previously [5] and is a simplification of a more complex model of receptor systems [24]. Time-dependent changes of states were defined as: ACh . We defined negative V shift as hyperpolarizing shifts. Negative V shift lead to an increased occupancy of the states U1 and B1, i.e. to a more hyperpolarized state of the receptor model. Likewise, positive V shift were defined as depolarizing shifts as they lead to increased occupancy of the states U2 and B2, i.e. to a more depolarized state of the receptor model. A total removal of voltage sensitivity was achieved for large shifts locking the M2R in its hyperpolarized state (V shift ≦ -150 mV) and its depolarized state (V shift ≧ +150 mV) (S3A Fig). A change in voltage affinity due to the presence of ACh was represented as a change in the transition rate. The transition rate coefficient α was multiplied with the constant K in the voltage dependent transition between B1 and B2. Similarly, the ACh affinity of the M2R was represented by the transition rate coefficients δ 0 and γ 0 .
Markovian Modeling of I K,ACh . A previously developed description for I K,ACh was used to calculate this current [23]. The description of I K,ACh includes its dependency on extracellular K + concentration [K + ] o and rectification: with the Nernst potential E K . We calculated the conductance of K ACh channels as the maximum conductance g K,ACh,max multiplied with the open probability O from the K ACh channel sub-model:

Measured data and model parameterization
The parameters of the model were determined by iterative stochastic optimization as previously described [5]. Model parameterization was implemented in Matlab R2017a (The Mathworks Inc., Natick, MA) and the Matlab Parallel Computing Toolbox. An error function based on the root mean squared differences of the measured versus simulated features of the activation protocol, deactivation protocol, and concentration response curve from [5] was minimized.
Measured features were based on whole-cell voltage-clamp experiments [5]. I K,ACh for the activation and deactivation clamp protocols was recorded in the presence of 0.1 μM ACh. The 3 best measured currents, see S1 Fig, from [5] of each protocol were then fitted to a monoexponential equation, averaged and then normalized to C off at -10 mV: with the constants A and C, and the time-constant τ of activation (on) and deactivation (off).
We choose this strategy, as extracting mono-exponential features and then averaging them introduces less error in the overall reconstructed behavior than averaging the signals themselves. The same fitting approach was used during the model parameterization.
Concentration response curves were measured as the resulting I K,ACh at either a holding potential of -100 mV or +50 mV at different ACh concentrations (Fig 5A and 5E). The currents were then normalized to the current measured at maximal ACh concentration for each holding potential.
A total of 8 simulated features f s,i , consisting of six features for the voltage clamp protocols (i.e. C on , C off , A on , A off , τ on , and τ off ) and two for the concentration response curves (i.e. I KACh, norm,50mV and I KACh,norm,-100mV ), were compared to their corresponding measured features f m,i (Fig 5B-5D and 5F-5H): Other components of the cost function were used to constrain the behavior of the model at specific ACh concentration and V m , and to ensure high open probability of the channel. Without bound ACh and V m of +50 mV, the state U2 was forced to be maximal. Respectively, U1 was forced to be maximal without bound ACh and a transmembrane voltage of -100 mV. Equivalently, with bound ACh, state B2 and B1 were forced to be maximal at -100 mV and 50 mV, respectively. Further, the sum of B1 and B2 was forced to be maximal with bound ACh.

Single cell simulations
For simulations of single cell electrophysiology, the new receptor-channel model was integrated in the publically available model of human SN cells [6]. The formulations for the effect of ACh on I CaL and I f were left unaltered throughout the experiments. The model was first exported to Matlab using OpenCOR (www.opencor.ws) and then modified accordingly. Numerical integration was performed by using the integrated ode15s formulation provided by Matlab. We measured BCL to characterize the rate of spontaneous activation of the simulated SN cell at varying concentration of ACh and V shift . The maximum conductivity g K,ACh,max was set 0.0022 μS to reproduce previously published heart rate slowing in the presence of 0−0.1 μM ACh [12]. Furthermore, we measured the APD at 90% repolarization (APD 90 ) to characterize the cardiac AP.

Model parameters
The stochastic parameterization yielded the model parameters (Table 3). In comparison to the parameterization of our previous model [5], the fit error of the features from the activation and deactivation protocols was reduced from 1.5 to 0.68 despite the additional error terms. Respectively the squared fit error for each feature is; C on : 0.06, C off : 0.025, A on : 0.11, Aoff: 0.014, t on : 0.04, t off : 0.18, I KACh,norm,50mV : 0.03, I KACh,norm,-100mV : 0.003. The total error of the features is equal to the square root of the sum of the respective squared fit errors. The corresponding modeled and measured current traces of the voltage protocols are shown in S2 Fig [5] currents, using the activation and deactivation voltage protocols from [5]. Artifacts in the measured traces, which were cropped for the monoexponential fitting, are shown in grey.