In the current paper we model the pulmonary vein AP using an extension of the 1998 Fenton-Karma model (Fenton and Cherry, 2008): the four variable Bueno-Orovio Cherry Fenton (BOCF) model. This is a monodomain phenomenological model of the human ventricular AP first introduced in Bueno-Orovio et al. (2008):

Here u represents the transmembrane voltage, J_fi, J_so and J_si are phenomenological summations of the fast inward, slow outward, and slow inward currents respectively. J_fi is effectively gated by the gating variable v, J_so is voltage gated, and J_si is effectively gated by the product of the gating variables w and s. D_BOCF is either a spatially dependant diffusion constant (under the assumption of isotropic diffusion), or a diffusion tensor (under the assumption of anisotropic diffusion). We always take initial conditions at the resting state, where [u(0), v(0), w(0), s(0)] = [0, 1, 1, 0]. A full description of this model can be found in Bueno-Orovio et al. (2008).

Given that the shape of the emergent electrical activity recorded on the lasso catheter may be constrained by the underlying structure and function of the AP, a propagating AP was simulated using the detailed biophysical Courtemanche model for the human atrium (Courtemanche et al., 1998; Imaniastuti et al., 2014; Labarthe et al., 2014) as a proxy for clinical AP data. A generic AP from the Courtemanche model was modified to account for the electrical remodeling associated with AF (Courtemanche et al., 1999) and used as the initial stimulus for the BOCF model with parameters as defined in the sample fitting code in the appendix of Bueno-Orovio et al. (2008). These parameter choices were then evolved using the Nelder-Mead Simplex Algorithm (Nelder and Mead, 1965) (implemented by MATLAB's fminsearch), by minimizing the root mean squared error between subsequent APs (see Figure 3). With a spatial resolution Δx = 0.2 mm, a diffusion constant of D_Court = 2.615 was necessary for the simulated wavefront to match the conduction velocity of 48 cm/s observed clinically (Labarthe, 2013). To eliminate any effects from boundary conditions or transients from the stimulus, the fit was performed at the point x = 10 mm on a tissue cable 20 mm long. A cycle length of 600 ms was used to match the clinical data. This process resulted in the parameter choices defined in Table 1.

Figure 3 shows the shape of the propagating APs under the above conditions using the BOCF model with parameters as in Table 1, alongside the Courtemanche AF model as described in Courtemanche et al. (1999). The important qualities reproduced were conduction velocity (indicated by the simultaneous spike), upstroke velocity, and AP duration.

For the case of anisotropic diffusion an asymmetric finite difference method was used to simulate the BOCF model (see van Es et al., 2014 for full details). Since in general the degree of anisotropy for an individual patient is unknown, we included the principal axes and eigenvalues of the diffusion tensor as additional parameters to be optimized by our fitting algorithm. Physiological studies place the anisotropy ratio between 2 and 10, (Koura et al., 2002; Xie and Zemlin, 2016), which were used as bounds in our algorithms. The initial principal axes were placed at 45 degrees to the x and y axes, maximizing the effect on the propagation pattern.

2D simulations of the pulmonary vein were performed by numerical integration of Equation (1) by a finite difference method over a discretized cylindrical domain to represent the excitable myocardial sleeve extending over the base of the pulmonary vein. Dimensions vary from vein to vein, with the right inferior typically the largest and the left inferior the smallest (Stojanovska and Cronin, 2008). We assume dimensions within the range of observed measurements: a length of 15 mm (Cronin et al., 2007) and a diameter of 12.5 mm (Cabrera et al., 2002; Kim et al., 2005). A spatial resolution of Δx = 0.2 mm was used to discretize this cylinder into a rectangular domain of 200 × 75 grid points. Periodic boundary conditions were used along the lines x = 1 and x = 200, whilst Neumann boundary conditions were used along the lines y = 1 and y = 75 (where x and y represent nodes on the grid). We set the conductivity to 0 to model the effect of lesions due to ablated tissue at the relevant points, following the approach introduced in Dang et al. (2005); Reumann et al. (2008); Tobon et al. (2010). As we are only concerned with the effect on the AP propagation from the ablation process, we do not require a model of the thermodynamic processes of the catheter itself (Berjano, 2006; Suárez et al., 2010).

A visual representation of this structure is shown in Figure 2, which shows the propagating intracellular potential with the lines y = 1, y = h_a, and y = h_r annotated (Figure 2C). An ectopic is initiated from a stimulus along the line y = 1; the edge of the myocardial sleeve furthest from the atrial junction. Virtual ablation is performed by introducing a line of lesions on the circle y = h_a such that D(x, h_a) = 0. Conduction gaps are modeled such that D(x, h_a) = D_BOCF (for conductive tissue on small segments of the circle y = h_a). Consequently, semi-circular wavefront(s) will form on the other side of the lesions. Although loosely based on the underlying mechanisms, the values of the obtained parameters are phenomenological, and fit to the available data to ensure an accurate simulation on the lasso catheter electrodes, not to provide an estimation of the real location of the conduction gap.

We simulate pulmonary vein recordings from the lasso catheter across n electrodes (where n is typically 10 or 20), on y = h_r, where h_r > h_a. The electrodes are assumed to be equally spaced d = 200/n apart, such that for an n electrode catheter c = (a, h_r) where a = {d, 2d, …, nd}). At each point c = (x′, y′), an approximation for the surface potential Φ described originally in Gima and Rudy (2002) is given by:

where

Bipolar recordings between electrodes i and j (denoted PV i-j clinically) are simulated by:

Throughout this paper, we divide the pulmonary vein into three equal sections, with the ablation line positioned at h_a = 25 and the recording catheter positioned at h_r = 50. This is a practical consideration, as quantifying these measurements during the procedure would be difficult given information collected as standard in clinical practice.

The important characteristics of both the simulated and recorded data are the activation times (from maximal absolute value of dV/dt) of each signal compared to the others, as this gives a representation of the wavefront shape termed the relative activation time curve. It is necessary to use the absolute value as the recordings are bipolar. The relative activation time curve can be visualized by plotting the catheter along the x axis and its activation time on the y axis, giving a curve of the activation times of each signal relative to the others.

To understand the relationship between the relative activation time curve and parameters of the overall pulmonary vein model, the quantity, size and locations of conduction gaps are used as input parameters, since these have the most profound effect on the emerging wavefront shape. The root mean squared error between relative activation times obtained from simulated and clinical recordings are minimized, again using the Nelder-Mead Simplex Search method (implemented by MATLAB's fminsearch) to establish the location of conduction gaps which result in the most accurate activation time curve. Here it is important to note we do not claim to have found the location of the conduction gap(s) via this fit, only that we have calibrated model parameters that most closely recreate the phenomenology of the waveforms from the recording catheter.

In the cases for which there is poor contact between recording catheter and tissue, the signal is typically flat or white noise. This is evident, for example, in channels PV 3-4 17-18, and 19-20 in Figure 1. To reconstruct missing electrode recordings, a partial relative activation time curve was obtained from the active channels. Model parameters of the overall pulmonary vein model were calibrated from the active channels, using the Nelder-Mead Simplex Search (implemented by MATLAB's fminsearch).

Pulmonary vein recordings used in this paper were obtained from adult male and female subjects undergoing pulmonary vein isolation therapy at Bristol Heart Institute. Bipolar recordings were obtained from a deflectable, circular, 20-pole Lasso catheter (Biosense Webster Ltd). Patients with both paroxysmal and persistent AF were included but all cases were paced into normal sinus rhythm by pacing at 600 ms intervals, as per standard clinical practice. All data were appropriately anonymized prior to their use in this study. Under United Kingdom law, patient data collected during normal clinical routine and anonymized before research use may be used for research without additional consent.

First we consider how well the model can reproduce the phenomenology of the pulmonary vein recordings when all 10 channels are active. To consider this, we use an exemplar set of clinical pulmonary vein recordings collected during pulmonary vein isolation therapy (as described in the methods). The goodness of fit between clinical recordings and model simulations is determined by minimizing the root mean squared error between the relative activation times of the model and the data. This is achieved by varying the positions of conduction gaps in the model. The average of ten recorded events in the data is used to form the target relative activation time curve. This ensures some robustness to variation in the data and enables us to estimate the conduction gap location and width as parameters, which should be constant until ablated.

For the chosen clinical data, and for parameter choices of the underlying BOCF model as in Table 1, we find that the root mean squared error between the relative activation time curve of the clinical data and that of the model (assuming isotropic diffusion) is minimized by placing conduction gaps centered on points x = 65.5 and x = 167, with widths 11 and 6 respectively. Both the number of minima and their locations are used to optimize position and width of the conduction gaps. This is important since both the number of minima and their locations within the relative activation time curve emerge as a result of the conduction gaps generating the signal. In current clinical practice, the earliest activation time(s) (e.g., the minima of the relative activation curve), are the most important, as these are assumed to be closest to the conduction gap and therefore the optimal ablation site. This is illustrated in Figure 4A, where we also present a model fit under the assumption of anisotropic diffusion. In this case the conduction gaps are centered on points x = 59.5 and x = 164.5, with width 11 and 13 respectively.

In Figure 4B, we present a comparison between the original choice of clinical pulmonary vein recordings and simulations for the two classes of model. Time units of the model are rescaled such that the relative activation time-scale of the model is equivalent to that of the clinical recordings, which permits a clearer visual comparison. Note that both classes of model result in visually similar simulated pulmonary vein recordings. We perform a more rigorous analysis of differences between anisotropic and isotropic diffusion later, when considering the ability of the model to reconstruct missing channels in the clinical data.

Next, we tested the capacity of the model to predict future ectopic events, given an average over an initial ten events. For the identified choice of model parameters from the initial ten events, we simulated a series of additional ectopic events and for each event we calculated the root mean squared error between relative activation time curves obtained from either simulated or clinical ectopic events. We define t₀ as the time of the last event in the training set, and t_N as the time of the N^th subsequent ectopic event. Figure 5 shows how this error scales as the number N of the ectopic event moves further away from the training set. The apparent periodicity in the error is most likely due to rhythmic movements of the patient, such as breathing.

We now focus on the capacity of the model to reconstruct missing channels, exemplars of which were shown in Figure 1. This is a key result of this paper, and the one with most relevance to a potential clinical decision support system. To test the accuracy of the model, we start with a clinical recording for which all channels are active. We then eliminated a subset n (n = 0 to 5) of the channels replacing them with a 0 time trace. Five was chosen as the upper limit, since clinically a recording with less than half the channels active would not be relied on for determining the site of ablation. We then estimated model parameters using the same approach as in the previous section, but only data from those channels that were active. Using these parameters we then contrasted the error between the relative activation time curves obtained from the simulated next ectopic event and the subsequent ectopic event from the original clinical recording (including all channels). This enables us to assess how well the model can reproduce clinically relevant information (since the relative activation time curve is used for determining the site of ablation).

Figure 6 shows a box plot for each value of n. Each box in the box plot represents the root mean squared error between the relative activation time curves obtained from the average across 20 model simulations (with anisotropic diffusion and without) and that obtained from a clinically recorded ectopic event. The case n = 0 enables us to consider the limit of the goodness of fit between the model and the clinical data. This is effectively the intrinsic error attributable to the choice of model. For subsequent plots, n random channels were removed from the training set (simulating the effect of lost information due to poor contact). Different time intervals and different signals were used for each calculation so that the error distribution presented is as close as possible to the errors that we might expect to observe in clinic. This is important as it minimizes the likelihood of observations simply being due to an artifact of the ectopic event chosen for the fit.

For up to 3 channels reconstructed, the median and maximum errors do not significantly increase over that of the control whether or not anisotropy is considered. This is an important result as it demonstrates reconstructing up to three missing channels is not a significant source of additional error and therefore the model as presented may ultimately have clinical use under these conditions. Removing more than 3 channels leads to information loss resulting in outliers with statistically significant errors (see the cases for 4 and 5 channels removed). This demonstrates the limit of the number of missing channels that the considered models can reliably reconstruct.

To consider the whether the assumption of anisotropic diffusion is significant, we performed an unequal variances t-test (so-called Welch's t-test) to test whether the errors from each model could have come from a distribution with the same mean. This test consistently showed no significant difference (p > 0.05), except for the case of 3 signals reconstructed (p = 0.0414). This suggests that whilst anisotropy is clearly important in terms of the underlying physiology, it does not significantly affect the quality of model fit to the phenomenology of the recorded signals. This is an important consideration as calibration of model parameters is more efficient under the assumption of isotropic diffusion.

To test the potential of this technique to aid the clinical procedure, we presented a cardiologist specializing in pulmonary vein isolation therapy, with three variations of clinical recordings collected from three patients:

the original clinical recordings with all channels active;

These scenarios are illustrated in Figure 7. The cardiologist was unaware of the origin of each recording, and to avoid bias, the recordings were supplied in a random order. The following results were obtained (summarized in Table 2 for convenience). For patient 1, given the original data, the first point of ablation would have been around PV 5-6, with PV 15-16 noted as a second choice. With channels 15-16, 17-18, and 19-20 removed, only PV 5-6 was identified as the only appropriate ablation zone. When these channels were reconstructed by the model, PV 17-18 was identified as the second choice of ablation target. For this patient, the model has helped to identify a second relevant ablation target that was not identified when channels were missing. If initial ablation is not successful, the cardiologist will ablate in the area surrounding the target area, hence an initial estimate closer to the optimal location will result in successful isolation using fewer radiofrequency pulses. This will result in a smaller region of tissue being ablated and a shorter procedure.

In the second patient presented, a clear earliest spike time was present on PV 5-6. The removal of PV 5-6 and its neighbors led to PV 9-10 being identified as an ablation target. In this case, the reconstruction led to the same zone being targeted as the original signal, while the estimation with the recordings missing was two channels away. For this patient, reconstruction of the missing recordings led the cardiologist directly to the optimal decision.

In the final patient, the earliest spike, on PV 7-8, was removed, along with PV 5-6 and 9-10. As previously, these missing channels shifted the chosen ablation target by 2 channels. The reconstructed signal led to a target selected which was closer to the target chosen with all information present. As with patient 1, we infer this result as satisfactory, as starting closer to the optimal target will lead to quicker isolation of the pulmonary vein.

In all three cases, the missing channels influenced the decision made by the cardiologist, demonstrating the potential impact of information loss in clinic. However, when the cardiologist used the recordings combined with the signals reconstructed by the mathematical model to make a decision, the decision made was closer to the decision that would have been made had all information been present. Whilst these results provide only limited proof of concept at this stage, assuming the original data and clinical decision to represent the “ground truth,” then we believe there is significant potential for our approach to minimize the effect of this lost information.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer NT and handling Editor declared their shared affiliation, and the handling Editor states that the process met the standards of a fair and objective review.