An atrioventricular node model incorporating autonomic tone

Abstract

The response to atrial fibrillation (AF) treatment is differing widely among patients, and a better understanding of the factors that contribute to these differences is needed. One important factor may be differences in the autonomic nervous system (ANS) activity. The atrioventricular (AV) node plays an important role during AF in modulating heart rate. To study the effect of the ANS-induced activity on the AV nodal function in AF, mathematical modelling is a valuable tool. In this study, we present an extended AV node model that incorporates changes in autonomic tone. The extension was guided by a distribution-based sensitivity analysis and incorporates the ANS-induced changes in the refractoriness and conduction delay. Simulated RR series from the extended model driven by atrial impulse series obtained from clinical tilt test data were qualitatively evaluated against clinical RR series in terms of heart rate, RR series variability and RR series irregularity. The changes to the RR series characteristics during head-down tilt were replicated by a 10% decrease in conduction delay, while the changes during head-up tilt were replicated by a 5% decrease in the refractory period and a 10% decrease in the conduction delay. We demonstrate that the model extension is needed to replicate ANS-induced changes during tilt, indicating that the changes in RR series characteristics could not be explained by changes in atrial activity alone.

1 Introduction

Atrial fibrillation (AF) is the most common supraventricular tachyarrhythmia (Hindricks et al., 2020). Characteristic for AF is an increased and irregular atrial activity that results in a rapid and irregular ventricular activation. Atrial fibrillation is linked to substantial morbidity and mortality, and is a significant burden to patients, physicians, and healthcare systems globally. Two main strategies of AF treatments are rate control and rhythm control. Rate control is one of the corner stones of AF management, however the effect of individual rate-control drugs are difficult to predict in advance. This is why the choice of a rate-control drug today remains empiric and driven largely by their safety profile and contraindications rather than predicted efficacy. Therefore, the complex mechanisms of AF have to be better understood to personalize the treatment and reduce the burden of AF on the healthcare system.

It has been shown that the autonomic nervous system (ANS) is contributing to the initiation and maintenance of AF (Shen and Zipes, 2014). Either a predominance in sympathetic or in parasympathetic modulation has been observed to initiate an episode of paroxysmal atrial fibrillation (PAF); and in some patients, both the sympatho-vagal and vagal predominances have been observed to initiate PAF episodes (Lombardi et al., 2004). Hence, differences in the ANS activity among patients may be an important factor behind the inter-patient differences in response to treatment. To investigate the ANS-induced changes to the pathophysiology of AF, the effect of the ANS has to be quantified. One common method to quantify the autonomic tone during normal sinus rhythm (NSR) is by heart rate variability (HRV) (Sassi et al., 2015). In sinus rhythm, HRV can be used to obtain information about the function of the sinoatrial (SA) node. This information is valuable for the quantification of the autonomic tone, because the SA node is densely innervated by the ANS (Shen and Zipes, 2014; George et al., 2017). In AF, however, HRV cannot be used to quantify the autonomic tone, because the heart beats are not initiated in the SA node.

Instead, the ventricular rhythm during AF is determined by the atrial electrical activity and the subsequent AV nodal modulation. Since the AV node is densely innervated by the ANS, characterizing the AV nodal behavior during AF may give valuable information about the autonomic tone. Results from previous studies suggest that the heart rate, as well as the heart rate variability, quantified by RR rmssd, and heart rate irregularity, quantified by RR sample entropy, are affected by β-blocker induced changes in sympathetic response (Corino et al., 2015). We hypothesize that such changes in the heart rate and its variability and irregularity reflect ANS-induced changes in the AV node. The ANS-induced changes on the cardiac electrophysiology can be studied using head-up and head-down tilt test, which in a previous study was shown to affect electrophysiological properties of atrial myocardium during AF (Östenson et al., 2017). It is unclear if the changes in the heart rate and its variability and irregularity are explained by the changes in the atrial electrophysiology alone or also by changes in the AV nodal properties. Investigating how the ANS is modulating the heart rate during AF is a complex task and requires a model based analysis.

Previously, several AV node models have been proposed that incorporate important characteristics of the AV nodal structure and electrophysiology in their design. Characteristic for the AV node is its dual-pathway physiology enabling a parallel excitation propagation of impulses with different electrophysiological properties (George et al., 2017). For example, the slow pathway (SP) has a longer conduction delay and shorter refractory period compared to the fast pathway (FP) (George et al., 2017). Furthermore, the refractory period and conduction delay are dynamic and depend on the recent history of the conducted and blocked impulses in the AV nodal tissue (George et al., 2017; Billette and Tadros, 2019). Early models of the AV node did not account for the dual-pathway physiology (Cohen et al., 1983; Jørgensen et al., 2002; Rashidi and Khodarahmi, 2005; Mangin et al., 2005; Lian et al., 2006). Later models have incorporated this feature, represented by separate refractory periods (Corino et al., 2011; Henriksson et al., 2016; Inada et al., 2017; Wallman and Sandberg, 2018) and separate conduction delays (Climent et al., 2011b; Inada et al., 2017; Wallman and Sandberg, 2018). However, no models have explicitly incorporated ANS-induced changes in their model description.

Therefore, the aim of the present study is to incorporate ANS-induced changes into the AV node network model previously proposed by Wallman and Sandberg, (2018). The extension of the AV node model was guided by a distribution-based sensitivity analysis (Pianosi and Wagener 2018) and incorporates ANS-induced changes in the computation of the refractoriness and conduction delay. The extended model is evaluated with respect to its ability to replicate changes in heart rate and RR series variability and irregularity observed during head-up and head-down tilt test.

2 Materials and methods

First, the clinical tilt test data is described in Section 2.1. The RR series characteristics are defined in Section 2.2, followed by the description of a network model of the AV node (Section 2.3). A sensitivity analysis on the AV node model is described in Section 2.4, that identifies the influence of changes in model parameters on the RR series characteristics. Based on the sensitivity analysis, the AV node model is modified to account for ANS-induced changes in AV node characteristics (Section 2.5). The ability of the modified AV node model to replicate ANS induced changes in RR series characteristics observed during tilt-test is assessed in Section 2.6. Finally, the statistical analysis is described in Section 2.7, that is used to determine significant differences in AFR and RR series characteristics between tilt positions.

2.1 Tilt test study

The autonomic influence on the RR series characteristics was analysed using ECG data recorded during a tilt test study performed by Östenson et al. (2017). Recordings from 24 patients with persistent AF were considered of sufficient quality for analysis and were included in the present study; their age was 66 ± 9 (mean ± std), and 63% were men. None of the patients had abnormal levels of thyroid hormones, severe renal failure requiring dialysis, or heart valve disease. None of the patients were ablated for AF or on any of the Class I or Class III antiarrhythmic drugs. The tilt test was performed between 1 and 3 p.m. in a quiet study room. Standard 12-lead ECG was recorded during supine position, followed by head-down tilt (HDT, -30°) and then head-up tilt (HUT, +60°). The tilt table was manually operated and had hand grip and ankle support for HDT and foot board support for HUT; the patients remained in each position approximately 5 min. ECG preprocessing and R-peak detection was performed using the CardioLund ECG parser (www.cardiolund.com).

2.2 RR series characteristics

The RR series consists of the intervals between consecutive heartbeats, where the time of a heartbeat is determined by the corresponding R peak in the ECG signal. In this work, three statistical measures of the RR series characteristics were used, quantifying heart rate, heart rate variability and heart rate irregularity, respectively, defined according to Eqs. 1–3. The mean of the RR intervals is computed as

where RR_i denotes the i:th RR interval in the RR series. The root mean square of successive RR interval differences (RR_V, variability) is computed as

The sample entropy of the RR series (RR_I, irregularity) is computed as

where the binary variable with , has the value 1 if the maximum absolute distance between corresponding scalar elements in the vectors , RR_i+1, …, and is below the tolerance r times the standard deviation of the RR interval series, otherwise the value is zero (Richman and Moorman, 2000). In this study, the parameters were set to m = 2 and r = 0.2.

2.3 Network model of the human atrioventricular node

The AV node is modelled by a network of 21 nodes (cf. Figure 1) (Wallman and Sandberg, 2018; Karlsson et al., 2021). The AV nodal dual-pathway physiology with a slow pathway (SP) and a fast pathway (FP) is represented with two chains of 10 nodes each. The last nodes of the two pathways are connected to each other and to an additional coupling node (CN). Impulses enter the AV node model simultaneously at the first node of each pathway and leave the model over the CN. Retrograde conduction is possible due to the bidirectional conduction within the pathways and between the last nodes of SP and FP.

FIGURE 1

Each node represents a section of the AV node and is described with an individual refractory period R^P(Δt_k) and conduction delay D^P(Δt_k) defined as

where P ∈ {SP, FP, CN} denotes the association to a pathway. The electrical excitation propagation through the AV node is modelled as a series of impulses that can either be passed on or blocked by a node. This decision is based on the interval Δt_k between the k:th impulse arrival time t_k and the end of the (k–1):th refractory period computed as

If Δt_k is positive, the impulse is conducted to all adjacent nodes, otherwise the impulse is blocked due to the ongoing refractory period R^P(Δt_k−1). The conduction delay D^P(Δt_k) describes the time delay between the arrival of an impulse at a node and its transmission to all adjacent nodes. If an impulse is conducted, R^P(Δt_k) and D^P(Δt_k) of the current node are updated according to Eqs. 4–6. For the computation of R^P(Δt_k) and D^P(Δt_k), the nodes in each pathway are characterized by six parameters, defining minimum refractory period, ; maximum prolongation of refractory period, ΔR^P; time constant ; minimum conduction delay, ; maximum prolongation of conduction delay, ΔD^P; and the time constant . The SP, FP and CN are modelled with separate vectors , ΔR^P, , , ΔD^P, ], all with fixed values.

The AV node model processes the impulse propagation chronologically and node by node, using a priority queue of nodes, sorted by impulse arrival time; details can be found in Wallman and Sandberg (2018). The input to the AV node model is a series of atrial impulses that is used to initialize the priority queue. As the impulses are conducted to adjacent nodes, new entries are added to the priority queue. The output of the AV node model is a series of impulses activating the ventricles.

In this study, the series of atrial impulses during AF is modelled as a point-process with independent inter-arrival times according to a Pearson Type IV distribution (Climent et al., 2011a). Hence, the atrial activation (AA) series is completely characterized by four parameters, namely the mean μ, standard deviation σ, skewness γ and kurtosis κ.

2.4 Distribution-based sensitivity analysis

The sensitivity of the three RR series characteristics to the AV node model and AA series parameters is evaluated by applying a distribution-based sensitivity analysis, based on the work of Pianosi and Wagener, (2018). For the sensitivity analysis, cumulative distribution functions (CDF) are estimated using a dataset of K = 250 000 randomly generated model parameter sets x and the characteristics y of the corresponding simulated RR series. For each simulation, an atrial impulse series with 60 000 AA intervals was generated using the Pearson Type IV distribution, with μ randomly drawn from  ms, σ randomly drawn from  ms, and γ and κ kept fixed to 1 and 6, respectively. The γ and κ were kept fixed since they cannot be estimated from the f-waves of the ECG. Negative AA intervals were excluded from the impulse series. The model parameters θ^SP and θ^FP were randomly drawn from bounded uniform distributions given in Table 1, as previously done in Karlsson et al. (2021). The θ^CN were kept fixed according to Table 1, corresponding to R^P(Δt_k) and D^P(Δt_k) of the CN equal to 250 ms and 0 ms, respectively.

The RR series characteristics were computed using a series of 4000 RR intervals corresponding to the first impulses that left the AV node model through the CN. Two selection criteria were used to remove non-physiological parameter sets. First, a model parameter set was only included if the slow pathway had a lower refractory period R^SP(Δt_k) < R^FP(Δt_k) and higher conduction delay D^SP(Δt_k) > D^FP(Δt_k) than the fast pathway for all Δt_k. Second, the resulting was required to be in the range 300  ms  1000 ms, corresponding to heart rates between 60 bpm and 200 bpm. Heart rates below 60 bpm are disregarded, because the pacemaker function of the AV node, that becomes relevant in this case (George et al., 2017), is not incorporated in the AV node model. Heart rates above 200 bpm are disregarded based on a reported minimum refractory period in the bundle branches of around 300 ms (Denes et al., 1974).

A sensitivity coefficient is computed for each pair of model parameter x_n and RR series characteristic y_m, where x_n is the n:th element in x and y_m is the m-th element in y. The indicates how much a change in model parameter x_n affects the distribution of y_m and is defined as

where is the Kolmogorov-Smirnov (KS) distance between the unconditional CDF and the conditional CDF . When estimating , the range of variation of x_n is split into C = 15 equally spaced conditioning intervals , with c = 1, …, C (cf. Figure 2A). All samples within are used to estimate the corresponding (cf. Figure 2B). To generate the set of , with d = 1, …, D, a subset of K/C samples are bootstrapped D = 1000 times (cf. Figures 2A,B). The KS distance is defined as

FIGURE 2

As the and are approximations based on a finite number of samples, parameters that have no influence on y_m can result in above zero. The impact of approximation errors on can be estimated for each y_m using a dummy parameter defined as

A model parameter x_n is determined to have influence on y_m if and only if .

2.5 Extended atrioventricular node model accounting for autonomic nervous system induced changes

The results from the sensitivity analysis (Section 3.1) indicate that changes in both the AV node model parameters and the AA series parameters have an influence on the RR series characteristics. Based on this, the AV node model described in Section 2.3 is extended to account for ANS-induced changes in the AA series by allowing μ(t) and σ(t) of the Pearson Type IV distribution to vary over time. Moreover, the AV node model was extended by two scaling factors A_R and A_D, accounting for the effect of changes in autonomic tone on refractory period (A_R) and on conduction delay (A_D).

The factors A_R and A_D model the combined effect of changes in sympathetic and parasympathetic activity and do not differ between the SP, FP and CN.

2.6 Tilt-induced changes in extended atrioventricular node model

In this section, the extended AV node model proposed in Section 2.5 is investigated with respect to its ability to mimic tilt-induced changes in RR series characteristics.

The clinical ECG signals (cf. Section 2.1) are used to generate AA series for the AV node model input and to compare the characteristics of the simulated RR series to the clinical RR series (cf. Figure 3). For this purpose, a continuous 15-min ECG signal with 5 minutes per supine, HDT and HUT position was desired for each patient. In the clinical data, however, the length of the three tilt positions varied between patients with the supine position being between 5 and 13 min, HDT being between 5 and 7 min and HUT being between 5 and 9 min. For two patients, there was an additional minute in supine position between the HDT and HUT. The ECG signals were aligned to the middle of the HDT section and a 15-min long segment centered around the same midpoint was chosen for each patient (cf. Figure 4).

FIGURE 3

FIGURE 4

The clinical RR series characteristic trends , and for each patient p are computed from the RR intervals using a sliding window of length N according to Eqs. 1–3 (cf. Figure 3). For RR_V and RR_I, N is set to 200, because shorter RR interval series might lead to inaccuracies in the sample entropy computation (Yentes et al., 2013). For , N is set to 100, as its computation is more robust than the computation of RR_V and RR_I and shorter RR interval series allow for a better temporal resolution. RR intervals in the clinical RR series preceding and following ectopic beats were excluded. For the computation of RR_I according to Eq. 3, vectors with excluded RR intervals were omitted. The RR series characteristic trends of each patient , and were averaged over all 24 patients to obtain population-averaged clinical trends , and (cf. Figure 3).

For the generation of the AA series, first, an atrial fibrillatory rate (AFR) trend is estimated from each 15-min ECG segment (cf. Figure 3). The AFR is estimated by fitting a complex sinusoidal model to the f-waves of the ECG, following spatiotemporal QRST cancellation, as described in Henriksson et al. (2018). From each of the resulting AFR trends, the AA series parameters μ(t, p) and σ(t, p) are estimated by the mean and standard deviation of 1/AFR using 1-min sliding windows; the resolution of the AFR trend is 0.02 s (cf. Figure 3). Then, μ(t, p) and σ(t, p) are averaged over all 24 patients, resulting in the population-averaged trends μ(t) and σ(t) (cf. Figure 3). Finally, the AA series is iteratively generated (cf. Figure 3). The first AA interval is drawn from the Pearson Type IV distribution with μ(0) and σ(0), and each consecutive AA interval is drawn from the distribution with μ(t_i) and σ(t_i) where t_i corresponds to the accumulated time of the previous AA intervals. The γ and κ of the Pearson Type IV distribution were kept fixed to 1 and 6, respectively.

For the simulations using the original and extended model, a set of 240 AV node model parameter vectors were generated (cf. Figure 3). Ten parameter vectors per patient were selected from a set of randomly drawn parameter sets based on their ability to replicate the RR series characteristics of the 5-min long supine segment of the respective patient. A detailed description of the parameter sets and the selection process can be found in the Supplementary Section 1. The ranges of the model parameters in the 240 parameter sets are given in Table 2.

For the computation of simulated RR series characteristic trends using the original and the extended model, respectively, simulations were performed with each of the 240 parameter sets using 10 different realizations of the AA series generated from μ(t) and σ(t). In the original model, the scaling factors A_R and A_D are not included, which is equivalent to the extended model using A_R = 1 and A_D = 1 (cf. Figure 3). In the extended model, A_R and A_D were allowed to change between supine and HDT and between HDT and HUT, respectively, but were assumed to remain constant within each position. Hence, for the extended model, A_R and A_D were set to 1 in the supine position, and different combinations of A_R ∈ {0.95, 1, 1.05} and A_D ∈ {0.8, 1, 1.2} were used for the simulations during HDT and HUT. For each simulation s, the mean RR interval trends and were computed from the RR interval series using a sliding window of length N = 100 (cf. Eq. 1). Whereas the RR variability and RR irregularity trends and , as well as and were computed from the RR interval series using a sliding window of length N = 200 (cf. Eqs. 2 and 3). The simulated RR series characteristic trends were averaged over all parameter sets and realizations to obtain the population-averaged simulated trends , and for the original model and , and for the extended model (cf. Figure 3).

2.7 Statistical analysis

A Wilcoxon signed rank test was applied to determine if AFR, , RR_V and RR_I differed significantly between supine, HDT and HUT. For the analysis, the AFR and RR series characteristics were computed for each patient and tilt position using the 5-min long ECG segments (cf. Figure 4). A p-value was considered significant.

3 Results

3.1 Sensitivity analysis

Results from the distribution-based sensitivity analysis (described in Section 2.4) with respect to the influence of the AV node model parameters on RR series characteristics are shown in Figure 5. Heart rate, quantified by is predominantly sensitive to changes in the refractory period parameters with the four largest contributors being the R_min and ΔR parameters of both pathways. In contrast, the changes in the conduction delay had little influence on the , with ΔD^SP being the only conduction delay parameter that is slightly above the dummy threshold. Changes in the mean of the AA series μ were also influential on the , while changes in the standard deviation σ of the AA series are not considered to have influence to changes in .

FIGURE 5

For RR_V quantifying RR series variability, nearly all model parameters of the refractory period, conduction delay and AA series had sensitivity coefficients above the dummy threshold. The four largest contributors to changes in the RR_V were the ΔR parameters of both pathways, as well as the minimum refractory period and minimum conduction delay of the fast pathway, and .

The RR_I quantifying RR series irregularity was also influenced by most model parameters of the refractory period, conduction delay and AA series. The four largest contributors were the minimum refractory period of the slow pathway , the standard deviation σ of the AA series and the maximum prolongation of the refractory period and conduction delay of the slow pathway, ΔR^SP and ΔD^SP.

3.2 Clinical data

The AFR decreased significantly from the supine position to HDT and increased significantly from HDT to HUT, where the AFR during HUT was significantly higher than during supine (Table 3). The heart rate increased during HDT and increased further during HUT (Table 3). The results align with the observations of Östenson et al. (2017). The variability and irregularity of the RR series decreased during HDT and decreased further during HUT (Table 3). For the variability and irregularity of the RR series, only the differences between supine and HUT were statistically significant.

3.3 Tilt-induced changes in atrioventricular node model

The average of μ(t) and σ(t) over all 24 patients is illustrated in Figure 6. The μ(t) shows a clear variation during HDT and HUT, but not in supine position, where μ(t) was approximately constant around 150 ms. Compared to μ(t) during supine position, μ(t) increased during HDT and decreased during HUT (cf. Figure 6).

FIGURE 6

In Figure 7, the characteristics , and estimated from clinical data during tilt test are illustrated. It can be seen that , and are decreasing from supine to HDT and decreasing further from HDT to HUT. When performing the simulations with the original model, , and are decreasing from supine to HDT, but increasing from HDT to HUT. When performing the simulations with the extended model, , and are decreasing from supine to HDT and decreasing further from HDT to HUT. Comparing the clinical and simulated trends of and RR_V(t), it can be seen that the extended model accounting for ANS-induced changes can better replicate the observed changes to the clinical RR series characteristics compared to the original model. For RR_I(t), both the original and extended model produce RR series that are more regular than the clinical RR series, as the irregularity quantified by the sample entropy is higher for the clinical RR series. For the simulated RR series characteristics of the extended model, the average of , and during the 5 min in HDT and HUT were illustrated for the nine different combinations of A_R ∈ {0.95, 1, 1.05} and A_D ∈ {0.9, 1, 1.1}. For , and , an increase in A_R causes an increase, and for and , an increase in A_D causes an increase. However, for , an increase in A_D instead causes a decrease. The , and are obtained using A_R = 1 and A_D = 0.9 during HDT and A_R = 0.95 and A_D = 0.9 during HUT and are displayed in Figure 7; the scaling factors were chosen so that the resulting and matches and .

FIGURE 7

4 Discussion

The aim of this study was to extend the AV node model (Wallman and Sandberg, 2018) to incorporate ANS-induced changes. The extension of the AV node model was guided by a distribution-based sensitivity analysis. The sensitivity analysis indicated that the refractory period and conduction delay parameters as well as the atrial impulse series had a significant influence on the heart rate as well as the variability and the irregularity of the RR series, while the most influential parameters were predominantly those describing the refractory period. Rather than modelling the effect of the sympathetic and parasympathetic activity separately, we describe the joint effects, i.e., the autonomic tone. We proposed an extension to the AV node model that accounts for the ANS-induced changes by introducing scaling factors for the refractory period and conduction delay. The capability of the extended AV node model to replicate ANS-induced changes was investigated by comparison to ECG data acquired during tilt test.

Our results (Figure 7) indicate that the extended model, but not the original, could replicate the observed changes in the clinical RR series characteristics during HUT and HDT, since the changes in RR series characteristics could not be explained by changes in atrial activity alone. The , and (Figure 7) show that a decrease in refractory period and conduction delay allow the model to replicate the decrease in , and . Conversely, if the refractory period and conduction delay are kept fixed for , and , all three RR series characteristics increase during HUT, which is the opposite direction of change of , and . When comparing and with , it can be seen that the sample entropy of the simulated RR series is lower than that of the clinical RR series. This highlights that the simulated RR series are more regular than the clinical RR series. One possible explanation for a lower irregularity in simulated RR series is the lack of short-term variations in AV node refractoriness and conduction delay. Such short-term variations may be induced by respiratory modulation in ANS activity. Thus, a natural next step in our model development will be to incorporate the respiratory modulation of the ANS, likely via periodical variations in the scaling factors A_R and A_D.

Many electrophysiological (EP) studies have demonstrated that an increase in sympathetic activity is causing a decrease in the human AV nodal conduction delay (Lister et al., 1965; Dhingra et al., 1973; Morady et al., 1988; Cossú et al., 1997) and a decrease in the refractory period (Morady et al., 1988; Cossú et al., 1997). Moreover, a decrease in sympathetic activity in the human AV node is causing an increase in conduction delay and refractory period (Morady et al., 1988). Head-up tilt is associated with increased sympathetic tone, and it has been demonstrated that the AV nodal conduction delay and refractory period decrease when changing the posture from supine to standing (Hashimoto et al., 1991). The results in Figure 7 confirm that a reduction in the conduction delay using A_D = 0.9 and a reduction in the refractory period using A_R = 0.95 better replicate the observed changes in the clinical RR series characteristics than the original model during HUT. Decreases in refractory period and conduction delay of up to 30% in response to isoproterenol-induced increases in sympathetic activity have been reported (Lister et al., 1965; Dhingra et al., 1973; Cossú et al., 1997). However, when considering that the reported changes in heart rate due to the isoproterenol administration is larger than the observed changes in during tilt, the parameter choice of A_D = 0.9 and A_R = 0.95 are reasonable for the tilt test data used in this study.

Increased parasympathetic activity has been associated with an increased conduction delay (Martin, 1977); studies in dogs reported an increased conduction delay with acetylcholine administration (Priola et al., 1983; Bertrix et al., 1984) and vagal stimulation (Spear and Moore, 1973; Martin, 1975; Pirola and Potter, 1990). Moreover, there are indications that an increased parasympathetic activity is associated with an increased refractory period (Martin, 1977); experimental studies using rabbit hearts reported an increased AV-nodal refractory period (West and Toda, 1967) and occurrences of 2:1 AV nodal block (Cranefield et al., 1959) with acetylcholine administration, and studies in dogs reported occurrences of AV block with acetylcholine administration (Hageman et al., 1985) and vagal stimulation (Spear and Moore, 1973; Hageman et al., 1985).

It is unclear how the HDT affects the sympathetic and parasympathetic activity. The results in Figure 7 show that a reduction in the conduction delay using A_D = 0.9 and no modification of the refractory period using A_R = 1 better replicate the observed changes in the clinical RR series characteristics than the original model during HDT. These results are consistent with possible slight increase in sympathetic tone provoked by HDT. However, other interpretations are possible. Nagaya et al. (1995) postulated a diminished sympathetic activity in HDT. Under that hypothesis, the results in Figure 7 suggest a decrease in parasympathetic tone to revert the direction of change caused by a decreased sympathetic tone. It should be noted that the model presented here does not distinguish between these two possibilities, since A_D and A_R are modelling the joint effect of changes in parasympathetic and sympathetic activity. Hence, the scale factor A_D = 0.9 during HDT could be reflecting either a slight increase in sympathetic activity, a slight decrease in parasympathetic activity, a larger increase in sympathetic activity combined with an increase in parasympathetic activity, or a large decrease in parasympathetic activity combined with a decrease in sympathetic activity.

The set of scaling factors A_R and A_D used to create , and in Figure 7 results in RR series characteristics similar to that observed during HDT and HUT. The results in Figure 7 show that a scaling factor A_R below 1, i.e., a decrease of the refractory period, causes a decrease in , and . Conversely, a scaling factor A_R above 1, i.e., an increase of the refractory period, causes an increase in , and . A scaling factor A_D below 1, i.e., a decrease in conduction delay, causes a decrease in , and vice versa. The opposite relationship can be seen for the RR series irregularity, where a scaling factor A_D below 1 causes an increase in and vice versa. Moreover, when considering one RR series characteristic at a time, it can be anticipated in Figure 7 that the same 5-min average value of the RR series characteristics can be achieved with different combinations of A_R and A_D. Hence, considering all three RR series characteristics simultaneously increases the likelihood of identifying a unique pair of scaling factors A_R and A_D that fits the observed data.

To reduce the complexity of the model, the refractoriness and conduction delay of the SP, FP and CN are modified with the same A_R and A_D. However, due to the structural and molecular heterogeneity of the different pathways, it is likely that the ANS-induced changes affect each pathway differently (George et al., 2017). In rabbit hearts, it was reported that acetylcholine strongly affects fibers of the atrionodal junction but does not show any effect in the lower part of the node or the bundle of His (Trautwein, 1963). In the description of the AV node model, the CN is merging the impulses from the SP and FP and its refractory period and conduction delay is independent of Δt_k. In contrast to Karlsson et al. (2021), was set to the minimum of the bounded uniform distributions for the and the given in Table 1. Further, the conduction delay of the CN was set to 0, as other choices of a constant conduction delay would not have changed the resulting RR series. In previous work on the network model (Wallman and Sandberg, 2018; Karlsson et al., 2021), the AA interval series was modelled as a Poisson process. However, based on results of Climent et al. (2011a), a Pearson Type IV distribution better reproduces the statistical properties of the AA interval series during AF and was therefore chosen in the present study. The mean and standard deviation of the Pearson Type IV distribution were determined from the mean and standard deviation of the AFR. However, the skewness and kurtosis were fixed, as their sensitivity coefficients were uninfluential (data not shown) and since there is no straight-forward way to estimate these parameters from the f-waves of the ECG.

In the present study, the ability of the extended model to mimic tilt-induced changes was investigated using data from a previous study (Östenson et al., 2017), with tilt angles fixed to -30° in HDT and 60° in HUT, respectively. Different tilt angles of the tilt, i.e., different magnitude of the orthostatic stimulus, may affect the ANS response and hence the resulting RR series characteristics. Previous results from patients in normal sinus rhythm show that the sample entropy of the RR series was decreasing during HUT from 0° to 60° but remained roughly constant from 60° to 90° (Porta et al., 2007). Based on these results, we assume that the tilt angle of 60° is sufficiently large to induce changes in autonomic tone. Access to data from patients with AF during other tilt-inclinations could potentially be used to refine the model to take the degree on inclination into account. The tilt-induced changes in RR series irregularity observed in the present study are in line with the results in Patel et al. (2018), where a decrease in RR sample entropy in response to HUT in patients with AF was reported. The tilt-induced changes in RR series irregularity observed in the present study are also in line with the changes reported for patients in normal sinus rhythm during HUT (Porta et al., 2007). Results from previous studies suggest that the RR series irregularity during normal sinus rhythm increase in response to HDT (Porta et al., 2015), whereas a slight but not significant decrease was observed in the present study with patients in AF. However, it should be noted that origin of RR series variability and irregularity during AF differs from that during normal sinus rhythm and hence, the interpretation of the results with respect to autonomic tone may be different.

The effect of the ANS-induced activity was investigated with respect to its ability to mimic the population-averaged changes observed during tilt test. The RMSSD and sample entropy were used to quantify RR series variability and irregularity, respectively, since these statistical measures have been used in previous studies to assess changes in RR series characteristics during AF in response to drugs (Corino et al., 2015) and tilt-test (Patel et al., 2018). Population-averaged trends were chosen over the trends of individual patients to reduce the uncertainty in the estimation of the clinical RR_V and RR_I trends. The parameter sets used for the simulations in Section 2.6 were selected to be representative of the patients in the present study based on their ability to replicate RR series characteristics observed during supine position. However, it should be noted that fitting of the model to individual patients is outside the scope of the present study. Due to the short measurement duration of the clinical data, a robust estimation of individual model parameters is not to be expected with the present methodology (Karlsson et al., 2021). Longer measurements from more patients will allow model development and evaluation on a patient-specific basis, forming an attractive next step.

A distribution-based sensitivity analysis was chosen over a variance-based method, because the distributions of the simulated RR series characteristics are highly-skewed and multi-modal. Hence, variance alone cannot adequately represent the uncertainty (Pianosi and Wagener, 2018). Instead, a distribution-based method characterizes the uncertainty and sensitivity by investigating the entire distribution of the model outputs (Pianosi and Wagener, 2018). The results of the sensitivity analysis in Figure 5 indicate that is the only model parameter that is uninfluential, since the sensitivity coefficients for all three RR series characteristics are below the dummy threshold. One important outcome of the sensitivity analysis therefore is that the refractory period and conduction delay of the AV node as well as the atrial input are influencing the RR series characteristics. For simplicity, we are proposing a linear scaling of refractory period and conduction delay parameters, but it would be interesting to refine this model description in the light of additional clinical data. It should be noted that the sensitivity coefficients are quantifying sensitivity on a global scale, and that there may be large local variations. As a result, the extent of variation in , RR_V(t) and RR_I(t) for a set of scaling factors A_R and A_D depend on the model parameters. For example, in Figure 7, it is clear that the scaling factor A_D affects , while the sensitivity analysis (Figure 5) indicates that the influence of changes in conduction delay on is very limited on a global scale.

In the present study, the estimates in and were based on sliding windows of N = 200 RR intervals. The choice of N is a tradeoff between estimation accuracy and time resolution. The sample entropy estimation is expected to stabilize with greater N and a minimum of N ≥ 200 was recommended by Yentes et al. (2013). In the present study, N was chosen as short as possible in favour of time resolution to investigate the ANS-induced changes in the RR series characteristics during tilt. To accommodate the estimation uncertainty resulting from a small N, the simulated RR series characteristics trends were averaged over 10 repeated simulations for 240 different parameter sets. For the sensitivity analysis, N was chosen to be 4,000 in favour of estimation accuracy since the simulation was stationary.

While ANS modulation has been extensively studied during normal sinus rhythm (Porta et al., 2007; Porta et al., 2015; Sassi et al., 2015; Patel et al., 2018), no attempts have been made towards the estimation of ANS modulation during persistent AF. The present study is a first step towards developing a model of the AV node that will ultimately be used to quantify ANS modulation on a patient specific basis by fitting to RR interval series and information on atrial electrical activity obtained from clinical ECG recordings. The results (Figure 7) show that the proposed extended model of the AV node accounting for changes in autonomic tone can better replicate changes in RR series characteristics observed during tilt-test than the original model, implying that this is a viable approach to take. Further developments are needed to incorporate ANS modulation in the model and methodology for robust estimation of such modulation from clinical data.

5 Conclusion

We present an extended AV node model that incorporates ANS-induced changes. The extension was guided by a distribution-based sensitivity analysis showing that changes in refractoriness and conduction delay of the AV node as well as changes in atrial activity significantly influence the RR series characteristics. We demonstrate that the model extension is needed to replicate the changes in heart rate and RR series variability and irregularity observed during head-up and head-down tilt.

References

• 1

  BertrixL.BouzouitaK.LangJ.LakhalM.ChahQ. T.FauconG. (1984). Potentiation by hypokalemia of the effects of acetylcholine on the canine heart in situ. Naunyn. Schmiedeb. Arch. Pharmacol.326, 169–174. 10.1007/BF00517315

  • CrossRef
  • Google Scholar

• 2

  BilletteJ.TadrosR. (2019). An integrated overview of av node physiology. Pacing Clin. Electrophysiol.42, 805–820. 10.1111/pace.13734

  • CrossRef
  • Google Scholar

• 3

  ClimentA. M.AtienzaF.MilletJ.GuillemM. S. (2011a). Generation of realistic atrial to atrial interval series during atrial fibrillation. Med. Biol. Eng. Comput.49, 1261–1268. 10.1007/s11517-011-0823-2

  • CrossRef
  • Google Scholar

• 4

  ClimentA. M.GuillemM. S.ZhangY.MilletJ.MazgalevT. N. (2011b). Functional mathematical model of dual pathway av nodal conduction. Am. J. Physiol. Heart Circ. Physiol.300, 1393–1401. 10.1152/ajpheart.01175.2010

  • CrossRef
  • Google Scholar

• 5

  CohenR. J.BergerR. D.DushaneT. E. (1983). A quantitative model for the ventricular response during atrial fibrillation. IEEE Trans. Biomed. Eng.30, 769–781. 10.1109/TBME.1983.325077

  • CrossRef
  • Google Scholar

• 6

  CorinoV. D.SandbergF.MainardiL. T.SörnmoL. (2011). An atrioventricular node model for analysis of the ventricular response during atrial fibrillation. IEEE Trans. Biomed. Eng.58, 3386–3395. 10.1109/TBME.2011.2166262

  • CrossRef
  • Google Scholar

• 7

  CorinoV. D.UlimoenS. R.EngerS.MainardiL. T.TveitA.PlatonovP. G. (2015). Rate-control drugs affect variability and irregularity measures of RR intervals in patients with permanent atrial fibrillation. J. Cardiovasc. Electrophysiol.26, 137–141. 10.1111/jce.12580

  • CrossRef
  • Google Scholar

• 8

  CossúS. F.RothmanS. A.ChmielewskiI. L.HsiaH. H.VogelR. L.MillerJ. M.et al (1997). The effects of isoproterenol on the cardiac conduction system: Site-specific dose dependence. J. Cardiovasc. Electrophysiol.8, 847–853. 10.1111/j.1540-8167.1997.tb00845.x

  • CrossRef
  • Google Scholar

• 9

  CranefieldP. F.HoffmanB. F.Paes de CarvalhoA. (1959). Effects of acetylcholine on single fibers of the atrioventricular node. Circ. Res.7, 19–23. 10.1161/01.RES.7.1.19

  • CrossRef
  • Google Scholar

• 10

  DenesP.WuD.DhingraR.PietrasR. J.RosenK. M. (1974). The effects of cycle length on cardiac refractory periods in man. Circulation49, 32–41. 10.1161/01.CIR.49.1.32

  • CrossRef
  • Google Scholar

• 11

  DhingraR. C.WinslowE.PougetJ. M.RahimtoolaS. H.RosenK. M. (1973). The effect of isoproterenol on atrioventricular and intraventricular conduction. Am. J. Cardiol.32, 629–636. 10.1016/S0002-9149(73)80055-4

  • CrossRef
  • Google Scholar

• 12

  GeorgeS. A.FayeN. R.Murillo-BerliozA.LeeK. B.TrachiotisG. D.EfimovI. R. (2017). At the atrioventricular crossroads: Dual pathway electrophysiology in the atrioventricular node and its underlying heterogeneities. Arrhythm. Electrophysiol. Rev.6, 179–185. 10.15420/aer.2017.30.1

  • CrossRef
  • Google Scholar

• 13

  HagemanG. R.NeelyB. H.UrthalerF.JamesT. N. (1985). Negative chronotropic and parasympatholytic effects of alinidine on canine sinus node and AV junction. Am. J. Physiol.248, 324–330. 10.1152/ajpheart.1985.248.3.H324

  • CrossRef
  • Google Scholar

• 14

  HashimotoT.FukataniM.MoriM.HashibaK. (1991). Effects of standing on the induction of paroxysmal supraventricular tachycardia. J. Am. Coll. Cardiol.17, 690–695. 10.1016/S0735-1097(10)80185-8

  • CrossRef
  • Google Scholar

• 15

  HenrikssonM.CorinoV. D.SörnmoL.SandbergF. (2016). A statistical atrioventricular node model accounting for pathway switching during atrial fibrillation. IEEE Trans. Biomed. Eng.63, 1842–1849. 10.1109/TBME.2015.2503562

  • CrossRef
  • Google Scholar

• 16

  HenrikssonM.MarozasV.SandbergF.SörnmoL. (2018). PetrnasModel-based assessment of f-wave signal quality in patients with atrial fibrillation. IEEE Trans. Biomed. Eng.65, 2600–2611. 10.1109/TBME.2018.2810508

  • CrossRef
  • Google Scholar

• 17

  HindricksG.PotparaT.DagresN.ArbeloE.BaxJ. J.Blomström-LundqvistC.et al (2020). 2020 ESC Guidelines for the diagnosis and management of atrial fibrillation developed in collaboration with the European Association for Cardio-Thoracic Surgery (EACTS): The Task Force for the diagnosis and management of atrial fibrillation of the European Society of Cardiology (ESC) Developed with the special contribution of the European Heart Rhythm Association (EHRA) of the ESC.Eur. Heart J.42, 373–498. 10.1093/eurheartj/ehaa612

  • CrossRef
  • Google Scholar

• 18

  InadaS.ShibataN.IwataM.HaraguchiR.AshiharaT.IkedaT.et al (2017). Simulation of ventricular rate control during atrial fibrillation using ionic channel blockers. J. Arrhythm.33, 302–309. 10.1016/j.joa.2016.12.002

  • CrossRef
  • Google Scholar

• 19

  JørgensenP.SchäferC.GuerraP. G.TalajicM.NattelS.GlassL. (2002). A mathematical model of human atrioventricular nodal function incorporating concealed conduction. Bull. Math. Biol.64, 1083–1099. 10.1006/bulm.2002.0313

  • CrossRef
  • Google Scholar

• 20

  KarlssonM.SandbergF.UlimoenS. R.WallmanM. (2021). Non-invasive characterization of human AV-nodal conduction delay and refractory period during atrial fibrillation. Front. Physiol.12, 728955. 10.3389/fphys.2021.728955

  • CrossRef
  • Google Scholar

• 21

  LianJ.MüssigD.LangV. (2006). Computer modeling of ventricular rhythm during atrial fibrillation and ventricular pacing. IEEE Trans. Biomed. Eng.53, 1512–1520. 10.1109/TBME.2006.876627

  • CrossRef
  • Google Scholar

• 22

  ListerJ. W.SteinE.KosowskyB. D.LauS. H.DamatoA. N. (1965). Atrioventricular conduction in man: Effect of rate, exercise, isoproterenol and atropine on the p-r interval. Am. J. Cardiol.16, 516–523. 10.1016/0002-9149(65)90028-7

  • CrossRef
  • Google Scholar

• 23

  LombardiF.TarriconeD.TundoF.ColomboF.BellettiS.FiorentiniC. (2004). Autonomic nervous system and paroxysmal atrial fibrillation: A study based on the analysis of RR interval changes before, during and after paroxysmal atrial fibrillation. Eur. Heart J.25, 1242–1248. 10.1016/j.ehj.2004.05.016

  • CrossRef
  • Google Scholar

• 24

  ManginL.VinetA.PagéP.GlassL. (2005). Effects of antiarrhythmic drug therapy on atrioventricular nodal function during atrial fibrillation in humans. Europace7, S71–S82. 10.1016/j.eupc.2005.03.016

  • CrossRef
  • Google Scholar

• 25

  MartinP. (1975). Dynamic vagal control of atrial-ventricular condition: Theoretical and experimental studies. Ann. Biomed. Eng.3, 275–295. 10.1007/BF02390973

  • CrossRef
  • Google Scholar

• 26

  MartinP. (1977). The influence of the parasympathetic nervous system on atrioventricular conduction. Circ. Res.41, 593–599. 10.1161/01.RES.41.5.593

  • CrossRef
  • Google Scholar

• 27

  MoradyF.NelsonS. D.KouW. H.PratleyR.SchmaltzS.De BuitleirM.et al (1988). Electrophysiologic effects of epinephrine in humans. J. Am. Coll. Cardiol.11, 1235–1244. 10.1016/0735-1097(88)90287-2

  • CrossRef
  • Google Scholar

• 28

  NagayaK.WadaF.NakamitsuS.SagawaS.ShirakiK. (1995). Responses of the circulatory system and muscle sympathetic nerve activity to head-down tilt in humans. Am. J. Physiol.268, R1289–R1294. 10.1152/ajpregu.1995.268.5.R1289

  • CrossRef
  • Google Scholar

• 29

  ÖstensonS.CorinoV. D.CarlssonJ.PlatonovP. G. (2017). Autonomic influence on atrial fibrillatory process: Head-up and head-down tilting. Ann. Noninvasive Electrocardiol.22, e12405. 10.1111/anec.12405

  • CrossRef
  • Google Scholar

• 30

  PatelH. C.HaywardC.WardleA. J.MiddletonL.LyonA. R.Di MarioC.et al (2018). The effect of head-up tilt upon markers of heart rate variability in patients with atrial fibrillation. Ann. Noninvasive Electrocardiol.23, e12511. 10.1111/anec.12511

  • CrossRef
  • Google Scholar

• 31

  PianosiF.WagenerT. (2018). Distribution-based sensitivity analysis from a generic input-output sample. Environ. Model. Softw.108, 197–207. 10.1016/j.envsoft.2018.07.019

  • CrossRef
  • Google Scholar

• 32

  PirolaF. T.PotterE. K. (1990). Vagal action on atrioventricular conduction and its inhibition by sympathetic stimulation and neuropeptide y in anaesthetised dogs. J. Auton. Nerv. Syst.31, 1–12. 10.1016/0165-1838(90)90166-g

  • CrossRef
  • Google Scholar

• 33

  PortaA.Gnecchi-RusconeT.TobaldiniE.GuzzettiS.FurlanR.MontanoN. (2007). Progressive decrease of heart period variability entropy-based complexity during graded head-up tilt. J. Appl. Physiol.103, 1143–1149. 10.1152/japplphysiol.00293.2007

  • CrossRef
  • Google Scholar

• 34

  PortaA.FaesL.MarchiA.BariV.De MariaB.GuzzettiS.et al (2015). Disentangling cardiovascular control mechanisms during head-down tilt via joint transfer entropy and self-entropy decompositions. Front. Physiol.6, 301. 10.3389/fphys.2015.00301

  • CrossRef
  • Google Scholar

• 35

  PriolaD. V.CurtisM. B.AnagnostelisC.MartinezE. (1983). Altered nicotinic sensitivity of AV node in surgically denervated canine hearts. Am. J. Physiol.245, 27–32. 10.1152/ajpheart.1983.245.1.H27

  • CrossRef
  • Google Scholar

• 36

  RashidiA.KhodarahmiI. (2005). Nonlinear modeling of the atrioventricular node physiology in atrial fibrillation. J. Theor. Biol.232, 545–549. 10.1016/j.jtbi.2004.08.033

  • CrossRef
  • Google Scholar

• 37

  RichmanJ. S.MoormanJ. R. (2000). Physiological time-series analysis using approximate entropy and sample entropy. Am. J. Physiol. Heart Circ. Physiol.278, 2039–2049. 10.1152/ajpheart.2000.278.6.H2039

  • CrossRef
  • Google Scholar

• 38

  SassiR.CeruttiS.LombardiF.MalikM.HuikuriH. V.PengC.-K.et al (2015). Advances in heart rate variability signal analysis: Joint position statement by the e-cardiology ESC working group and the European heart rhythm association co-endorsed by the asia pacific heart rhythm society. Europace17, 1341–1353. 10.1093/europace/euv015

  • CrossRef
  • Google Scholar

• 39

  ShenM. J.ZipesD. P. (2014). Role of the autonomic nervous system in modulating cardiac arrhythmias. Circ. Res.114, 1004–1021. 10.1161/CIRCRESAHA.113.302549

  • CrossRef
  • Google Scholar

• 40

  SpearJ. F.MooreE. N. (1973). Influence of brief vagal and stellate nerve stimulation on pacemaker activity and conduction within the atrioventricular conduction system of the dog. Circ. Res.32, 27–41. 10.1161/01.RES.32.1.27

  • CrossRef
  • Google Scholar

• 41

  TrautweinW. (1963). Generation and conduction of impulses in the heart as affected by drugs. Pharmacol. Rev.15, 277–332.

  • Google Scholar

• 42

  WallmanM.SandbergF. (2018). Characterisation of human AV-nodal properties using a network model. Med. Biol. Eng. Comput.56, 247–259. 10.1007/s11517-017-1684-0

  • CrossRef
  • Google Scholar

• 43

  WestT. C.TodaN. (1967). Response of the A-V node of the rabbit to stimulation of intracardiac cholinergic nerves. Circ. Res.20, 18–31. 10.1161/01.res.20.1.18

  • CrossRef
  • Google Scholar

• 44

  YentesJ. M.HuntN.SchmidK. K.KaipustJ. P.McGrathD.StergiouN. (2013). The appropriate use of approximate entropy and sample entropy with short data sets. Ann. Biomed. Eng.41, 349–365. 10.1007/s10439-012-0668-3

  • CrossRef
  • Google Scholar