The model of the AV node, shown in Figure 1, consists of a network of nodes and is based on the model presented in Wallman and Sandberg (2018). The model consists of two pathways, representing the SP and the FP, connected with a coupling node. Each pathway is modeled with 10 nodes, where each node corresponds to a localized part of the AV node. Each node can block incoming impulses or send them through adding a conduction delay. All nodes but the coupling node sends impulses to all other nodes connected to it, whereas the coupling node only receives impulses. A new refractory period [R_i(n)] and conduction delay [D_i(n)] are calculated every time a node (i) receives a new impulse (n). The refractory period and conduction delay are based on the stimulation history of the node and are described using exponential functions. These exponential functions have previously been used to fit AV node characteristics (Shrier et al., 1987; Lian et al., 2006; Wallman and Sandberg, 2018), and can be seen in Equations (1–3).

Here t~i(n) refers to diastolic interval preceding impulse n, t_i(n) the arrival time of impulse n at node i, and t_i(n − 1) and R_i(n − 1) the arrival time and refractory period of impulse n − 1 at node i, respectively. If t~i(n) is negative, the node will still be in its refractory period and thus the impulse will be blocked. The model parameters defining minimum refractory period, R_min; maximum prolongation of refractory period, ΔR; time constant τ_R; minimum conduction delay, D_min; maximum prolongation of conduction delay, ΔD; and the time constant τ_D, are assumed to be fixed for the nodes in the SP and FP, respectively.

The coupling node models the total refractoriness and conduction delay introduced by the connection between the AV node and the bundle of His, the Purkinje fibers, and the bundle of His. This node has a separate set of parameters, representing separate functional properties, and will be denoted the His and Purkinje (HP) node. The refractory period for the Purkinje fibers is assumed to not affect the ventricular activation during AF. Thus, the whole refractory period for the HP node is determined by the bundle of His. However, the conduction delay for the HP node is viewed as the time it takes an impulse to travel from the start of the bundle of His to the end of the Purkinje fibers. The conduction delay from the start of the bundle of His until the end of the Purkinje fibers has clinically been showed to have a mean of 60 ms with a standard deviation of 10 ms for patients suffering from AF (Deshmukh et al., 2000). Thus, the conduction delay for the HP node is fixed at 60 ms. The HP node's refractory period is estimated by the mean of the ten shortest RR intervals, RR_min.

This results in 12 free parameters for the proposed model, denoted as a parameter vector θ=[RminFP, ΔRFP, τRFP, RminSP, ΔRSP, τRSP, DminFP, ΔDFP, τDFP, DminSP, ΔDSP, τDSP]. It is assumed that the first node of each pathway is simultaneously stimulated for incoming impulses from the atria. The model can then be used to produce a RR interval series with minimal computational demands using a modified version of Dijkstra's algorithm (Wallman and Sandberg, 2018). A link to the code for the model together with a basic user example can be found at section 5. The total minimum conduction delay and maximum prolongation, defined as Dmin,totFP=NnDminFP; ΔDtotFP=NnΔDFP; Dmin,totSP=NnDminSP; ΔDtotSP=NnΔDSP; where N_n = 10 are the number of nodes in each pathway, are introduced for convenience of presentation.

This study was based on ambulatory ECG data from the RATe control in Atrial Fibrillation (RATAF) study, which is approved by the regional ethics committee and the Norwegian medicines agency and was conducted in accordance with the Helsinki Declaration (Ulimoen et al., 2013). The RATAF study contains 24-h Holter recordings of 60 patients under baseline and during treatment with four different rate reducing drugs. All patients had permanent AF, no heart failure or symptomatic ischemic heart disease, an age of 71±9 (mean ± std), and 70% were men. To evaluate the presented model, we selected 15 min ECG segments, one for each of five patients, obtained under baseline conditions between 1:00 and 3:00 pm. These five patients were selected to be representative for the whole data set, with varying RR interval series characteristics and an average heart rate ranging between 63 and 140 bpm. In addition, corresponding ECG data obtained during treatment with Diltiazem was also used for one of the five patients.

The RR interval series were extracted from the ECG signals by first detecting the R peaks, before removing RR intervals preceding and following ectopic beats identified based on heartbeat morphology (Lagerholm et al., 2000). Along with this, the mean arrival rate of the atrium-to-atrium (AA) intervals was estimated from the f-waves in the ECG by first extracting the atrial activity from the ECG using spatiotemporal QRST cancellation (Stridh and Sornmo, 2001), before tracking the atrial fibrillatory rate (AFR) using a method based on a hidden Markov model (Sandberg et al., 2008). Finally, correction of the atrial fibrillatory rate by taking the atrial depolarization time into account was used to obtain an estimate of the arrival rate. Here, we denote the true mean arrival rate λ, and the estimated mean arrival rate λ^. One value of λ^ was obtained for each ECG segment (Corino et al., 2013).

Simulated data were created by fitting the model to the RR interval series from the five patients, cf. section 2.5, and using the resulting estimated model parameters to simulate an RR interval series of 20 min. The sequence of atrial impulses arriving to the AV node, and thus the input to the model, were simulated using a Poisson process with the mean arrival rate set to the value of λ^ estimated for each patient (Corino et al., 2011; Henriksson et al., 2015). The parameter values used for the simulated data, along with average heart rate of the simulated RR interval series, are summarized in Table 1.

To evaluate how well the model matches the extracted RR interval series, a fitness function comparing the model output to the RR interval series is used. In order to take the dynamics of the RR interval series into account, the Poincaré plot is used as a basis for the fitness function. The Poincaré plot is a scatter plot of successive pairs of RR intervals. To use the Poincaré plot as a fitness function, the RR interval series is binned into two dimensional bins centered between 250 and 1,800 ms in steps of 50 ms, resulting in N = 961 bins. The error function is computed according to Equation (4).

Here ϵ is the error value, and xi~ and x_i the number of RR intervals, in the i-th bin, of the measured data and model output, respectively. The normalization by x~i is introduced to avoid bins with a large number of data points to dominate the optimization. The square root is used as a trade-off between no normalization, making the bins with a large number of data points dominate, and normalization with the whole measured bin counts, making the accuracy of every bin have the same weight regardless of how much of the data are in that bin. A schematic representation of the parameter estimation process can be seen in Figure 2.

An initial study of how ϵ varies with varying model parameter values revealed a highly chaotic structure with a large number of local minima. This prompted us to minimize ϵ using a genetic algorithm (GA). A brief description of the algorithm is given below, with more detailed information in the Supplementary Section 1. Due to the high dimensional parameter space and the risk of premature convergence early in the optimization, a variant of an island model was used (Wahde, 2008). A schematic representation of the GA is shown in Figure 3. As visible in the figure, the full GA can be divided into two sections. The first section consists of five separate GA. This was implemented by restarting the algorithm five times with 300 individuals in each generation. The individuals in each starting run were initialized using a latin hypercube sampling in the ranges: {RminSP,RminFP}∈[250,600] ms; {ΔR^SP, ΔR^FP} ∈ [0, 600] ms, {τRSP,τRFP}∈[50,300] ms; {DminSP,DminFP}∈[0,30] ms; {ΔD^SP, ΔD^FP} ∈ [0, 75] ms; {τDSP,τDFP}∈[50,300] ms. These starting runs last for six generations, and after each run the best 150 of the individuals are saved and used in the second section, the main GA. Thus, the main GA uses a population of 750 individuals in each generation. For both the starting runs and the main GA, the 2.5% fittest individuals in each generation survives into the next generation unchanged, whereas the remaining individuals are created via tournament selection, two-point crossover, and creep mutation (Wahde, 2008). In order to avoid premature convergence, both incest prevention in the form of mating restriction between too similar individuals during crossover, and a varying mutation rate depending on the diversity of the individuals in each generation were implemented (Wahde, 2008). This process of selection, crossover, and mutation is then continued until termination. The termination of the starting runs always occurs after six generations. The termination for the main GA occurs either when ϵ for the fittest individual in each generation does not change for three generations, or when 15 generations have been run. The fittest individual for the k-th generation, ϵk^, is deemed to have changed if the difference between ϵ^¯(k) and ϵ^¯(k-2), seen in Equation (5), is lower than 25.

As described in section 2.3, a Poisson process with mean arrival rate λ^ was used as input to the model, and due to the stochastic nature of the Poisson process, ϵ varies between realizations. The magnitude of this variation was analyzed by finding a parameter set replicating the extracted RR interval series from patient 3 well, before simulating that parameter set with different lengths of the resulting RR interval series, L_RR, as seen in Figure 4, left panel. Each L_RR was simulated 1,000 times. Moreover, six more parameter sets with increasing ϵ were also simulated 1,000 times with the same L_RR, as seen in Figure 4, right panel.

The ϵ variation is decreasing with larger L_RR, however, the running time for the model is linearly increasing with L_RR, and thus shorter outputs are preferable. The variation of ϵ is not as important early in the optimization since the variation relative ϵ is smaller for larger ϵ, see Figure 4, right panel. However, after several generations most of the ϵ for individuals found by the GA are low, and thus the variability in ϵ has a larger impact on the algorithm. Therefore, L_RR is increased throughout the optimization.

As seen in Figure 3, the L_RR for all generations in the starting runs were 1,000 impulses. For the main GA, the first five generations used a L_RR of 3,000 impulses, the following five generations a L_RR of 5,000 impulses, followed by three generations with length of 7,500 impulses, before ending with two 10,000 impulses long generations. To obtain a robust estimate of ϵ^¯(k), the individual with the best fit in each generation is evaluated again with a L_RR of 10,000. After termination for the main GA, the 15 fittest individuals were tested again, with a L_RR of 50,000; this in order to select the fittest individual with a low variation in ϵ.

The effect of using an island based start together with varying L_RR was evaluated by comparing it to using only the main GA, as described in section 2.5, with L_RR fixed at 5,000. The initialization for this fixed GA was the same as for the starting runs, a latin hypercube sampling in the same ranges, and the population size was again 750. Performances of the two methods were evaluated by comparing the error value of the fittest individual for each generation, ϵ^k with the cumulative L_RR used for the evaluations, i.e., the accumulated total number of impulses in each generation. For the different starting runs, all runs were computed in parallel so that ϵ^k during this stage is the lowest value out of all the five starting runs. The average results from comparing the two versions of the GA on all five patients, each 100 times, are shown in Figure 5. From this it is possible to see that a lower ϵ^k, and thus a better fit to the RR interval series, can be found in less computational time using the proposed methodology. For reference, estimating the parameters for one patient using a single core on a standard desktop computer (Intel® Core™ i7-6600U Processor, @ 2.60GHz) requires on average 20 min, with variations due to the different terminating requirements for the GA. It is also possible to see that the termination criteria for a maximum number of generations stated in section 2.5 is typically achieved after the GA has converged.

Simulated RR interval series were used to evaluate the robustness of the model parameter estimates. The results from optimizing the model 200 times for the five simulated RR interval series can be seen in Table 2, where the mean and standard deviation for each of the 12 estimated parameters, for each of the five patients, are listed. Moreover, the mean error, defined as the difference between the mean value of the estimated parameter and the ground truth, averaged over the five patients, are also listed. Furthermore, the mean and standard deviation of the error normalized with respect to the parameter ranges, cf. section 2.5, are presented. From the SP ratio it is evident that the SP is used more for transmission, and from the normalized error, it is evident that the parameters associated with the SP are more robustly estimated. The histogram and Poincaré plots for the five simulated patients with the transmission pathway for each RR interval marked out can be seen in Supplementary Section 3, together with the simulated histograms showing the effect of changes to λ.

To set the robustness in perspective, the AV nodal properties were estimated 200 times for a single patient during baseline and under the influence of the non-dihydropyridine calcium channel blocker rate control drug Diltiazem. The results, shown in Figure 6, indicate that the uncertainty in the parameter estimation is sufficiently low in order to reveal the drug effect.

To evaluate the ability of the model and proposed workflow to represent AF data and to have a frame of reference, the proposed model is compared with the model presented in Wallman and Sandberg (2018); henceforth denoted the reference model. Both models were fitted to the RR interval series from one example patient, and the properties of the resulting simulated RR interval series are shown in the form of histograms, Poincaré plots, and autocorrelations, as seen in Figure 7. For both models, the optimizer was run until no change in error value for the fittest individual during ten generations occurred, to assure convergence. Both models used the optimizer described in section 2.5, but the reference model uses a fitness function based on the histogram (Wallman and Sandberg, 2018). It is clear from both the Poincaré plots and the autocorrelation plots that the proposed model can better replicate the dynamics of the RR interval series. The fit to the Poincaré plot can be quantified by the resulting ϵ, which for the proposed model was 1,360, compared to 6,740 for the reference model. Similarly, the value for the first lag autocorrelation was −0.07 for the proposed model and 0.52 for the reference, compared to the ground truth at −0.07.

The main limitation of the present study is the lack of comparison between the estimated parameter and the ground truth AV node characteristics, making the results more difficult to evaluate. Although simulated data was used as a substitute, it is not fully known how closely it matches reality. Another limitation is the assumption that both pathways are activated simultaneously, an assumption that may not be valid, since the electrical activity in the atria is highly disorganized. The variation in output originating from the stochastic input sequence can also be seen as a limitation to the proposed model, since the output for a single set of parameters can vary depending on the realization of the input sequence. However, without electrical measurements in the atria, it is not possible to model the exact behavior of the AV node.

Moreover, due to the computational time of estimating the parameters for each simulated RR interval series 200 times, only a subset of RATAF was used. However, the five patients were selected to ensure a representative subset based on their RR interval series characteristics. It should be noted that the focus of the present study is to evaluate the robustness in parameter estimation rather than analysis of the RATAF data set. Using the model to analyze the entire RATAF data set, including all patients, drugs, and time segments for outcome prediction forms a natural next step in this line of inquiry, and efforts toward this goal are ongoing at the time of writing.

Example results, cf. Figure 6, suggest that the estimates of refractory period and conduction delay are sufficiently robust to detect changes in response to treatment with rate control drugs. However, this needs to be verified in a larger study population. By using the model to simulate the treatment effect of different drugs in a patient-specific setting, it might be possible to predict the outcome of the drug treatment and thus assist in treatment selection. Furthermore, it could also be useful in drug development, by aiding in understanding what AV node properties are affected by a novel compound, and in what way.