Modeling specific action potentials in the human atria based on a minimal single-cell model

We present an effective method to model empirical action potentials of specific patients in the human atria based on the minimal model of Bueno-Orovio, Cherry and Fenton adapted to atrial electrophysiology. In this model, three ionic are currents introduced, where each of it is governed by a characteristic time scale. By applying a nonlinear optimization procedure, a best combination of the respective time scales is determined, which allows one to reproduce specific action potentials with a given amplitude, width and shape. Possible applications for supporting clinical diagnosis are pointed out.


Introduction
Detailed reaction-diffusion models to describe human atrial electrophysiology were first developed in the late 1990s [1][2][3][4] and are further developed until now. Important steps forward have been made to include specific ionic currents [5][6][7][8][9][10], which in particular allow one to investigate specific effects of pharmaceuticals in treatments of atrial fibrillation and other heart failures. Complementary to these detailed models, Bueno-Orovio, Cherry and Fenton introduced in 2008 a minimal reaction-diffusion model (BOCF model) for action potentials (AP) in ventricular electrophysiology, where the large number of ionic currents through cell membranes is reduced to three net currents [11]. This model has four state variables, one describing the transmembrane voltage (TMV), and the other three describing the gating of ionic currents. The TMV, as in detailed reaction models, satisfies a partial differential equation of diffusion type with the currents acting as source terms, and the time evolution of the gating variables is described by three ordinary differential equations coupled to the TMV. By fitting the action potential duration (APD), the effective refractory period and the conduction velocity to the detailed model of Courtemanche, Ramirez and Nattel [1] (CRN model), the BOCF model was recently adapted to atrial electrophysiology (BOCF model) [12].
In this work we develop a method to model specific AP based on the BOCF model as it is aimed in the clinical context in connection with improved and extended possibilities of diagnosis [13]. Compared to the detailed models, the BOCF model has the advantage that it is better amenable to some analytical treatment. This allows us to identify a small set of relevant model parameters for capturing the main features of a specific AP. Our methodology is a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 sketched in Fig 1 and can be summarized as follows. We start by labeling each given AP with its amplitude APA and with four APD, namely at 90%, 50%, 40% and 20% repolarization, denoted as APD 90 , APD 50 , APD 40 , and APD 20 respectively. These APD n (n = 20, 40, 50, 90) together with the amplitude APA are suitable to catch a typical shape of a specific AP, see Fig 2. The APD n taken for a specific patient are given to a parameter convertor that retrieves specific parameter values of the BOCF model. As relevant parameters, we adjust three time scales governing the closing and opening of the ionic channels. The parameter convertor consists of an optimization algorithm that searches for the best set of parameter values consistent with the measured AP properties.
The paper is organized as follows. In Section "BOCF model for atrial physiology" we shortly summarize the BOCF model and discuss the role of the three fit parameters that we selected to model specific AP. In Section "Parameter dependence of BOCF action potentials" we show how these parameters can be adjusted to obtain a a faithful representation of the AP properties APA, APD n , and in Section "Modeling of patient-specific action potentials with the BOCF model" we demonstrate the specific AP modeling for surrogate data generated with the CRN model [1]. A summary of our main findings and discussion of their relevance is given in Section "Conclusions". In the supporting information, we provide analytical calculations for the BOCF model that motivated our choice of fit parameters for the AP modeling. We also analyze the robustness of the optimization procedure with the activation frequency.

BOCF model for atrial physiology
The BOCF model has four state variables, which are the scaled TMV u, and three variables v, w and s describing the gating of (effective) net currents through the cell membrane. The TMV V is obtained from u via the linear relation V = V R (1 + αu), where for atrial tissue we set V R = −84.1 mV for the resting potential and α = 1.02 [12]. The time-evolution of u is given by the single-cell action potential model, here defined as where J = J(u, v, w, s) is the total ionic current and J stim an external stimulus current. The total ionic current decomposes into three net currents, a fast inward sodium current J fi = J fi (u, v), a slow inward calcium current J si (u, w, s), and a slow outward potassium current J so = J so (u), where the nonlinear functions F, G and H, are specified in S1 Appendix. There we show that the four differential Eqs (1) and (3) can be reduced to a system of two differential equations. This reduction shows that the three characteristic times τ fi , τ si and τ so1 , which fix the typical duration of the respective currents, govern the shape of the AP [cf. S1 Appendix]. We take these three time scales as parameters for fitting a specific AP and keep all other parameters fixed. For the values of the fixed parameters we here consider the set determined for the electrically remodeled tissue due to atrial fibrillation [12,14].

Parameter dependence of BOCF action potentials
In this section we show that in the BOCF model the amplitude APA can be expressed by a quadratic polynomial of the times τ fi , and the APD n by cubic polynomials of τ si and τ so1 . The dependence of APA and the APD n on the characteristic times, was determined from generated AP in single-cell simulations of the BOCF model by applying periodically, with a frequency f = 3 Hz, a square stimulus current of 40 pA, corresponding to an amplitude of 4.76 s −1 for the current J stim in Eq (1), for a time period of 3.5 ms. The resulting time evolution of the TMV in response to this stimulus was calculated by integrating Eqs (1) and (3)  , the APA depends only very weakly on τ si and τ so1 . Neglecting these weak dependencies, on τ si and τ so1 , we find the APA to increase monotonically with τ fi in the range [85, 110] mV relevant for human atria. In Fig 3(c) we show that the parameter τ fi can be well described by the quadratic polynomial where the coefficients c i and the coefficient of determination R 2 of the fit are given in Table 1. Likewise, as demonstrated in where the coefficients c ðnÞ mk are listed in Table 1 together with the R 2 values of the fits. The APD n of the single cell BOCF model depend on the activation frequency f or basic cycle length BCL = 1/f. Corresponding restitution curves are shown in Fig 5 for the remodelled tissue. These curves resemble the restitution curves known for others atrial models, see Ref. [15]. With higher frequency (shorter BCL) the APD n become smaller. This decrease is more pronounced for frequencies above 6 Hz. As a consequence, the optimization procedure becomes less robust for f ≳ 6 Hz, a feature that is discussed in more detail below in Section "Robustness with respect to the activation frequency".

Modeling of patient-specific action potentials with the BOCF model
Let us denote by V the APA and by D n the values of the APD n of a specific patient. To model the corresponding AP with the BOCF model, we determine τ fi by inserting APA ¼ V in Eq (4) and (τ si , τ so1 ) by minimizing the sum of the squared deviations between the the APD n , i. e. the function For the numerical procedure we used the Levenberg-Marquardt algorithm [16]. As one sees from Fig 4(b)-4(e), the APD vary monotonically with the time scales in the ranges fixed above. We checked that the Hessian is positive definite in the corresponding region, implying unique solutions when minimizing F . To demonstrate the adaptation procedure, we generated surrogate AP with the CRN model [1] for electrically remodeled tissue due to atrial fibrillation [14]. Specifically, we consider the maximal conductances, g Ca and g Na of the calcium and sodium currents to vary, while keeping all other parameters fixed to the values corresponding to the electrically remodeled tissue. The conductance g Ca affects both the AP plateau and the repolarization phase and the g Na controls mainly the amplitude of the AP [1].  from their values γ Na = 7.8 nS/pF and γ Ca = 0.0433 nS/pF for the electrically remodelled tissue [14]. The corresponding AP modeled with the BOCF model, i. e. for τ fi from Eq (4), and τ si and τ so1 obtained from the minimization of F ðt si ; t so1 Þ in Eq (6), are shown as dashed lines in the figures. In all cases these reproduce well the AP shapes generated with the CRN model.
To quantify the difference between the AP, we denote by A CRN ðtÞ and A BOCF ðtÞ their time course, and compute their relative deviation based on the L 2 -norm, where The initial time t i and final time t f are defined as the times for which u(t i ) = u(t f ) = θ 0 with θ 0 = 0.015473 (see S1 Appendix), with opposite signs of the corresponding time derivatives, i.e. du dt j t i > 0 and du dt j t f < 0. Fig 7(a) shows that, when keeping g Na = γ Na fixed, DA is below 5% for values of g Ca between 10-400% of the reference value γ Ca . For g Ca =g Ca ≳ 4, DA starts to increase. Likewise, as show in Fig 7(b), DA does not exceed 9% when varying g Na between 10-400% of γ Na for fixed g Ca = γ Ca .
Additionally to the relative deviation between AP, one can compute the relative deviations between the APA and APD n retrieved from the BOCF fit, Here X represents either V, giving Δ APA or D n , giving Δ APD n .  Fig 7(c) and 7(d) show ΔX APA as a function of g Ca /γ Ca and g Na /γ Na , again for fixed g Na = γ Na and g Ca = γ Ca , respectively. Corresponding plots of the ΔX APD n are shown in Fig 7(e) and 7(f) .  Fig 7(c) shows that ΔX APA is always very small, even for large deviations of g Ca from the reference value γ Ca . By contrast, ΔX APA is quite sensitive to variations of g Na . The deviation becomes larger than 5% for g Na =g Na ≳ 2.
As for the ΔX APD n they are typically below 12% except in the case of APD 20 . The APD 20 refers to the TMV level closest to the maximum and exhibits larger deviations up to about 20% for even small shape deviations.
All in all, Fig 7 shows that the optimization procedure retrieves acceptable fits of single-cell AP in a wide range of calcium and sodium conductances.  Table 1). The AP were taken after 10 4 ms, beyond the time needed for achieving the stationary state. https://doi.org/10.1371/journal.pone.0190448.g005

Robustness with respect to the activation frequency
The optimization procedure described in this paper was illustrated using one single activation frequency, namely f = 3 Hz. An important issue is the robustness of the optimization framework for other activation frequencies, which we address in Figs 8 and 9. Fig 8(a)-8(c) show the variation of the three coefficients in Eq (4) to fit the functional dependence of the APA on the parameter τ fi . As one sees, all three coefficients are approximately constant for activiations below 7 Hz. In that range of values one also observes a coefficient of determination R 2 ≳ 0:999, as shown in Fig 8(d). For higher frequencies f, the coefficients start to vary and the R 2 values of the fits become smaller, indicating the need of higher order polynomials to describe the relation between τ fi and APA. Similar results are obtained for the coefficients used to fit the APD surfaces as functions of the parameters τ si and τ so1 . These are shown in Fig 9 and demonstrate that the optmization procedure can be applied faithfully in the range 1-6 Hz. Outside this range, polynomials of higher order would be needed for better matches.
All in all, this section provides evidence that our optimization procedure derived for an activation frequency of 3 Hz, may also be applicable for frequencies ranging at least between 1 and 6 Hz. Fig 10 shows two illustrative examples of real AP and the respective fit with the optimization procedure. For details about the real data see Ref. [17].

Conclusions
In this work we showed how to model patient-specific action potentials by adjusting three characteristic time scales, which are associated with the net sodium, calcium and potassium ionic currents. The framework explores the possibilities of parameter adjustment of an atrial physiology model, namely the BOCF model [11], to reproduce AP shapes with a given amplitude, width and duration. The BOCF model is defined through a reaction-diffusion equation, coupled to three equations for gating variables that describe the opening and closing of ionic channels. It is simple enough to guarantee low computational costs for even extensive simulations of spatio-temporal dynamics [18]. Through a semi-analytical approach given in the S1 Appendix we showed why the three ionic currents suffice to derive the main features of empirical AP.
The high flexibility for case-specific applications can be used for clinical purposes. By adjusting a simulation to specific patient conditions one may also analyze numerically the effect of drug therapy under specific conditions. Using the optimization procedure for AP shape adjustment, the three characteristic times are retrieved, which are directly connected to the ion-type specific net currents. AP shapes showing pathological features will be reflected in Modelling case-specific action potentials the values of one (or more) times outside acceptable ranges. Accordingly, one can associate a corresponding net current and therefore identify the class of membrane currents, where pathologies should be present. In this sense the clinical diagnosis can be supported by the modeling. A future application could be to take the retrieved parameters values as a basis for spatially extended simulation by including the diffusion term in Eq (1) [11]. For this, one would need access to conduction properties which then would enable one to model spatiotemporal AP evoluion.
Though our framework is applicable in a quite wide range of values of sodium and calcium conductances, for conductances beyond a few times the reference values for electrically remodelled tissue the matching of AP shapes becomes less accurate. As for changes of the activation frequency f, the analysis in the Supporting Information limits the applicability of the AP modeling based on Eqs (4) and (5) to the range f = 1-6 Hz.
Furthermore, in case information is obtained about AP shapes from different places of the atria, e. g. by using a lasso catheter, a corresponding AP shape modeling would allow one to Modelling case-specific action potentials construct a patient-specific model with spatial heterogeneities. Based on this, it could become possible to generate spatio-temporal activation pattern and to identify possible pathologies associated in the dynamics of the action potential propagation. ranges of u-values, where the evolution of the set of variables changes discontinuously. Vertical dashed lines intersect the AP at one specific dotted line, thus bounding the time intervals corresponding to each region of u-values. In several of such time intervals, some of the variables decay exponentially and independently from the other variables, which simplifies the model considerably. In the regions where no exponential evolution is indicated the model follows the reduced system of equations derived in S1 Appendix. (EPS) S1 Appendix. Dynamical features of the BOCF model: A semi-analytical approach.

Supporting information
(TEX)