Biophysical models of cardiac electrophysiology (5) are an important tool to understand how cardiac tissue properties affect the generation and propagation of cardiac electrical signals (action potentials, APs). They also offer an ideal means for the training and development of computational tools that may aim to infer EP properties from electrical and optical mapping measurements, such as EP-PINNs.

Several biophysical EP models have been proposed (see models.cellml.org/electrophysiology), each with varying degrees of detail aiming to reproduce different EP features, cardiac regions or animal/human experimental findings. Mathematically, these EP models usually take the form of a reaction-diffusion system where a diffusion term or equivalent (6) models the propagation of the electrical signal across the cardiac tissue. In the monodomain formulation, a partial differential equation (PDE) describes the spatio-temporal variations in the electrical potential across a myocyte cell membrane (V). This PDE is usually coupled to one or more ordinary differential equations (ODEs) describing how, at each point in time and space, V and other local state variables both determine and are determined by the flux of ions across the cell membrane (5).

The most parsimonious model of the action potential describes it as a travelling excitation wave followed by a non-excitable (refractory) region. This representation requires at least two state variables: V, which spreads (diffuses) across neighbouring regions, and a non-observable, non-diffusible recovery variable W which effectively controls the refractoriness and restitution properties of the model. One of the simplest models that captures these properties is the 6-parameter canine ventricular Aliev-Panfilov model (7), which models the transmembrane ionic currents (V − W relationship) using smooth, differentiable functions. Furthermore, the diffusion of V across the cardiac tissue can be described by the monodomain equation (5), which, when combined with the Aliev-Panfilov model gives:

The diffusion term ∇→.(D∇→V) reduces to D∇²V in the case of homogeneous and isotropic conduction, i.e., when the diffusion tensor D is approximated by the same scalar throughout. Intuitively, this term quantifies how fast V is able to spread to its immediate neighbourhood to become more spatially homogeneous. D is determined mostly by the electrical conductivity of the myocardium and is a strong determinant of the propagation velocity of the AP. To prevent a leakage of V to regions outside the heart domain, the system described by Equations (1, 2) usually obeys no flux Neumann boundary conditions: ∂V∂n→=0 in the boundary of the heart tissue.

The −kV(V − a)(V − 1) − VW term in Equations (1, 2) models the rapid changes in V caused by ionic fluxes across the cell membrane. a is related to the excitation threshold (i.e., the minimum V value that leads to the onset of an AP). The model's APD and refractoriness can, in turn, be controlled using b. The values for each of the model parameters are typically chosen empirically to reproduce observed electrical signals - we use the values listed in Supplementary Table 1. The Aliev-Panfilov model uses rescaled units: V is adimensional (typically in the [0, 1]AU interval) and time is measured in temporal units, referred to as TU throughout this study. 1 TU corresponds to approximately 13 ms (7).

Other more complex EP models exist, through which it is possible to model individual membrane ionic currents and other biological components relevant for the AP and its propagation. One example used in the current study is a 14-current 30-variable canine atrial model that incorporates different degrees of EP remodelling caused by atrial fibrillation (8).

By assigning different sets of parameters and/or initial conditions to EP mathematical models, they can represent the electrical behaviour of the heart in both healthy and arrhythmic conditions. Healthy conditions are usually represented as unidirectional smooth propagation from a single source (Figures 1a,b). Large sources produce wave fronts which are close to planar, whereas point-like wavefronts lead to centrifugal (convex) wavefronts. Arrhythmias are usually modelled as one or more re-entrant waves (called spiral waves or rotors in 2D–Figure 1c). Moreover, the wavefronts and/or wavebacks of spiral waves can, in some instances, fragment (break-up), leading to complex activation patterns (9). These models can also consider localised pathology such as fibrosis, scar or ischaemia as heterogeneities in one or more model parameters. Localised reductions in D, for instance, can be used to reproduce the slow-down of AP propagation in fibrotic lesions (10). These heterogeneities can lead to local changes in the curvature of the activation wavefront (Figures 1d,e).

Biophysical models are usually employed in forward mode, with the aim of reproducing the system's behaviour assuming that all model parameters are perfectly known. In many circumstances, such as the identification of pathology from EP measurements, it is more desirable to use these biophysical models in inverse mode, inferring the tissue parameters that underlie an observed system behaviour. This task, often called parameter estimation (or systems identification), is extremely challenging for several reasons.

The observed data are typically insufficient to identify a unique parameter value, since the observations are usually sparse, incomplete and polluted by noise. This can be handled by optimisation methods, such as least squares, that fit data and model in an optimal way, combined with problem-dependent regularisation terms to stabilise the estimate (11). Typically, the optimisation requires many forward runs of the underlying model. Unfortunately, in physical systems that are described by PDEs such as in EP, the forward runs are very computationally expensive in realistic 2D and 3D settings. These difficulties are exacerbated by the fact that the parameters in many EP models are heterogeneous. Hence, we often do not infer a scalar quantity but a (discrete) function in space and/or time. This increases the cost and complexity of the inverse mode even further.

Modern inverse problem solvers such as ensemble Kalman filters (12) sequential Monte Carlo (13) and parameter estimation based on Markov chain Monte Carlo (14) are often combined with reduced order models (15) or multi-fidelity approaches (16) to decrease the computational complexity of parameter estimation. These very complex methods typically make strong assumptions about the statistical distribution of model parameters and require dedicated problem-specific parameterisation. It is not clear how well they can generalise when applied to a different EP model or task.

Physics Informed Neural Networks (PINNs) (17) are an exciting new tool for the study of physical systems modelled by PDEs and/or ODEs. PINNs have been shown to both efficiently find high-resolution solutions (forward problem) and perform parameter estimation (inverse problem) in a variety of systems (18). As opposed to most types of neural networks (NNs), whose inputs are exclusively empirical data, PINNs incorporate explicit knowledge about the physical laws that govern a system. This allows PINNs to compute solutions to initial and/or boundary value problems with comparatively less training data than conventional NNs.

Very briefly, PINNs make use of automatic differentiation (19) to rewrite the differential equations that a system obeys as the minimisation of a functional f. For example, for Equations (1, 2), this functional would be defined as:

PINNs are trained to minimise a hybrid loss function, L, which ensures the system obeys known physical laws (described by f_V and f_W), whilst simultaneously fitting known empirical system measurements. L thus includes terms to account for:

Mathematically:

Each of the terms of the loss function is typically computed in different domains:

There is an asymmetry in L between measurable (V) and latent (W) variables: only experimental measurements (and initial conditions) for V are usually available. Moreover, as W does not diffuse across the tissue, it does not obey any boundary conditions.

PINNs have typically been used for two main purposes (18). In the so-called forward mode, the NN's parameters are optimised to provide a representation of the physical system of interest consistent with observations. Using this representation, the system's differential equations can subsequently be solved at an arbitrarily high spatial or temporal resolution, bypassing the constraints (e.g., small temporal and spatial steps) of traditional numerical solvers. In inverse mode, PINNs additionally perform parameter inference (systems identification) by having the NN optimise one or more of the equation parameters (which here represent tissue EP properties) during the training process.

PINNs have recently been used in several areas of cardiovascular medicine, especially for applications related to blood flow. Examples include the estimation of myocardial perfusion and related physiological parameters from dynamic contrast enhanced MRI (20) and the estimation of haemodynamic parameters from microscopic images of aneurysms-on-a-chip (21).

In the field of cardiac EP, Sahli-Costabal et al. (22) used PINNs in forward mode to estimate activation time maps (ATs, i.e., the arrival times of the action potential) and conduction velocity (CV) maps in the left atrium at high spatial resolution. Sahli-Costabal's method uses PINNs to solve the (isotropic diffusion) eikonal equation, a simple relationship between ATs and the spatial gradient of CV. This effectively interpolates AT and CV across the left atrial surface. Although their PINNs implementation was exclusively deployed on simulated data, the proposed application is very clinically relevant, as it aims to mitigate the low spatial resolution of clinical AT measurements. Grandits et al. (23) subsequently extended this PINNs-eikonal equation framework to anisotropic conduction, using it to estimate high-resolution AT maps and fibre directions from in silico and patient data. The PINNs method performance was nevertheless lower than that of a traditional (variational) inverse solver (24). As they rely on the eikonal equation, these tools are not well suited to the study of arrhythmic conditions or to the inference of EP parameters other than AT and CV.

Maps of transmural electrical potential (V), similar to those simulated using the Aliev-Panfilov model, can be experimentally recorded using optical mapping. Optical mapping is a technique in which voltage-sensitive fluorescent dyes are added to cardiomyocyte preparations before imaging at high spatio-temporal resolution (25). It can be used to effectively obtain uncalibrated measurements of V(x→,t) in cardiac tissue across time. Although optical mapping can be challenging in vivo, in vitro experiments can provide very detailed insights into AP properties and cellular-level EP properties and gain insights into arrhythmic mechanisms (26). In particular, optical mapping can be used to study the effect of anti-arrhythmic drugs on cardiomyocyte preparations (27). These data were used to test EP-PINNs in an experimental setting.

We use an FD solver to generate in silico cardiac EP data, V_GT, which we use to train and evaluate the performance of EP-PINNs. Simulations were carried out in two different geometries: in 1D, in a 2-cm 1D domain cable; and in 2D, in a square with a 1-cm side. We use in-house code written in Matlab 2020b (Mathworks, Natick, MA, USA) that relies on central FD and an explicit 4-stage Runge-Kutta method to solve the isotropic monodomain Aliev-Panfilov model, as defined in Equations (1, 2) and model parameters listed in Supplementary Table 1. All simulations use Neumann boundary conditions and set dt = 5 × 10⁻³TU and dx = 100μm as temporal and spatial steps. All simulations were run for 300 TU, with the calculated VGT(x→,t) field saved at every 1 TU and at every spatial step (every 100μm). We thus generate in total 1.4 × 10⁴ and 7.0 × 10⁵ data points for VGT(x→,t), respectively in 1D and 2D.

Figure 1 shows example time frames of the generated V maps in 2D. All APs were initialised by adding an external stimulus current V_stim = 0.12 AU to the right-hand side of Equation (1) for 1 TU in a sub-domain of the studied geometry (Figure 1). In 2D simulations, this includes both planar (Figure 1a) waves (emanating from a rectangular stimulus) and centrifugal waves (from a point-like excitation) (Figure 1b). Spiral waves were also created using the cross-field protocol: as an initial planar excitation propagates, a second planar excitation wave, orthogonal to the first, is initiated. When timed appropriately (42 TU after the first stimulus, in our model), the second planar wave continuously curves as it moves to non-refractory tissue, giving rise to a sustained spiral wave (Figure 1c).

The data generated by the FD model is treated as GT data and is used to both train and test EP-PINNs. As in past studies (18, 28), we train EP-PINNs on a case-by-case basis, with a small subset of the data for which parameter inference is going to be performed. This is in contrast to most supervised NNs, which are usually trained with large amounts of data acquired in varied circumstances.

Training data for EP-PINNs, VGT(x→,t), are provided to the network as an input and used to minimise the data discrepancy loss function, L_data in Equation (6). They amount to 10–20% of the total generated data, corresponding to a variable number of data points as detailed below and in Supplementary Table 2. The remaining data are withheld from EP-PINNs and used, in a post-processing step, to assess EP-PINNs' ability to reproduce cardiac APs, as detailed below.

EP-PINNs are designed, trained and deployed using the Python DeepXDE library (28). As in past successful implementations of PINNs (18, 28), we use a fully connected network architecture.

As detailed in Figure 2, EP-PINNs take as inputs the spatio-temporal points, (x→,t), where they will estimate the main outputs: V (and W). Experimental measurements of V, VGT(x→,t), are also provided to EP-PINNs as inputs in a (training) subset of (x→,t). EP-PINNs minimise the hybrid physics-informed loss function described before (see Equations 5 and 6), by adjusting the network's weights and biases (collectively named θ in Figure 2). In inverse mode, EP-PINNs additionally adjust one or more parameters of the Aliev-Panfilov model (generically λ in Figure 2). The number of layers and neurons used by EP-PINNs is adjusted to the domain size and the complexity of the problem at hand, as described below and detailed in Supplementary Table 2.

The optimisation approach used for EP-PINNs also depends on problem size. For 1D problems, we use Adam optimisation (29)—see Figure 2C. We empirically determined that, in 2D, EP-PINNs' performance improved when initially using Adam optimisation for only the data agreement term (L_data in Equation 6), followed by Adam optimisation for the full loss function and ending in a final phase of L-BFGS optimisation (30)—see Figure 2C. The initial Adam training phases are used to rapidly approach the desired minimum and the final L-BFGS optimisation phase helps the network converge faster towards it (28).

In the presence of spatially-varying EP parameters (D in the current study), we use network architecture B (see Figure 2B). Here, D(x→,t) is estimated by a parallel NN, NN_D, with the same number of layers and neurons as the main NN. In this setup, D is treated as a system variable (on par with V and W) and its estimates directly contribute to the loss term that ensures the agreement with the EP equations: L_fV in Equation (6).

The hyperbolic tangent function (tanh) is used throughout as the differentiable activation function and Glorot initialisation from a uniform distribution is used for all weights (31). Additionally, to minimise convergence problems caused by explosive gradients (28) and enhance the NN's stability, we implemented an automatic reset of the training process when the losses at the first epoch of the training exceeded a predefined threshold.

EP-PINNs were trained on a high performance machine with 1 RTX6000 GPU and 4 AMD EPYC 7742 CPUs. Typical training times varied between 15 min (for 1D problems) and 16h (for heterogeneous spiral wave problems), as detailed in Supplementary Table 2.

We withheld the value of one or more of the model parameters from EP-PINNs, which were instead estimated by it. These parameters were chosen for their comparatively simple biophysical interpretation, known susceptibility to both disease remodelling and pharmacological action, and the limited degree of mathematical coupling between them. They were:

We tested EP-PINNs performance on in vitro datasets using optical mapping data from neonatal ventricular rat myocyte preparations stained with a voltage-sensitive dye, as described in detail by Chowdhury et al. (27). Briefly, we used four time series (movies) of optical mapping images (field of view: 4.1 × 0.1 mm, spatial resolution: 1.2 μm, temporal resolution: 2 ms, duration: 300 ms). In two of these image series, ionic channel modulating drugs (E-4031 or nifedipine) had been administered at half maximal inhibitory concentration (IC₅₀): 772.2 nM for nifedipine and 243.4 nM for E-4031. The other two temporal image series consisted of matched control (baseline) preparations, to which no drug had been given.

We manually selected two square regions of interest (ROIs) with a side of 2.3 μm and with their centres 0.7 mm apart, in the same image location for each time series. We spatially averaged the optical signal over each ROI to obtain a signal trace across time. From this signal, we manually selected two consecutive APs and normalised the signal to the [0, 1] interval for consistency with the Aliev-Panfilov model. To improve the signal to noise ratio (SNR) of this trace, we applied a mean average filter twice, aligned and averaged the two APs over time to obtain a single higher-SNR AP. These pairs of post-processed APs were used as inputs to 1D EP-PINNs in inverse mode, with b as the variable to be estimated. All model parameters were unchanged from those in Supplementary Table 1. We used EP-PINNs' architecture A with training scheme 2 to estimate b 10 times for each setting. 116-200 points were used for training EP-PINNs and 29–50 to test it, as detailed in Supplementary Table 2.

We investigated in particular whether EP-PINNs could detect the effect on APD of E-4031 and nifedipine, which, respectively, block the hERG voltage-gated potassium channel (I_Kr current) and the L-type calcium channel (I_CaL current). Whereas, E-4031 will extend APD (and thus shorten b in the Aliev-Panfilov model—see Equations 1, 2), nifedipine will have the opposite effect, decreasing APD (increasing b). We compared EP-PINNs' b estimates for i) E-4031 vs. baseline and ii) nifedipine vs. baseline and assessed, using a t-test, whether the administered drugs significantly changed b values.

Using optical mapping signals, EP-PINNs were able to accurately reproduce experimental APs and identify the actions of nifedipine and E-4031, correctly estimating that they, respectively, reduce and increase APD (Figure 9). The reduction in APD caused by nifedipine, an I_CaL blocker, was detected by EP-PINNs as a significant increase in b in the Aliev-Panfilov model [Δb = (2.7 ± 0.9) × 10⁻², p < 10⁻⁵]. The E-4031-induced increase in APD was more subtle [Δb = (−2.1 ± 3.0) × 10⁻²] and non-significant (p = 0.38). The reduced effect of E-4031 in these data is consistent with the modest role I_Kr, the current blocked by E-4031, is expected to play in rodent APs (32).

EP-PINNs offer a flexible and easy to implement framework for parameter estimation in EP. Sahli-Costabal et al. (22) and Grandits et al. (23) had already demonstrated PINNs' potential in cardiac EP by estimating high-resolution left atrial AT and CV maps in sinus rhythm conditions using a simple activation-only biophysical model. We extend PINNs' applications in EP by applying them to a more complex biophysical model, the monodomain Aliev-Panfilov model (7), which also captures restitution properties through the inclusion of a latent (non-measurable) variable, W. For the first time, we use a PINNs framework for the simulation of arrhythmic conditions, such as spiral waves, and for the estimation of parameters unrelated to the AT of V.

Importantly, we show EP-PINNs' are able to reproduce APs and perform parameter estimation in the absence of any data for W, which is not available experimentally. This bodes well for the deployment of PINNs for even more complex EP models, which use a higher number of latent variables to model individual ionic channels. This possibility is also supported by the work of Yazdani et al. (33), who successfully used PINNs for parameter estimation across several biological systems described by large sets of coupled ODEs.

We demonstrate PINNs' ability to describe AP dynamics in several circumstances. In 1D, EP-PINNs were able to reproduce APs even in the presence of very reduced amounts of experimental data (Figure 3A) and large amounts of noise (Figure 3B). PINNs' incorporation of explicit biophysical equations in the NN's loss function acts as an effective regulariser in EP problems, as demonstrated before in many other physical systems (17, 18). We note that these inherent regularisation properties allow PINNs to be trained with much lower amounts of training data than conventional NNs. The main drawback is the need for a comparatively time-intensive training on a case-by-case basis, compared to the global training usually employed with supervised NNs.

In 2D, we were able to accurately replicate AP dynamics for planar, centrifugal and planar waves, even in the presence of heterogeneities in the diffusion coefficient (see Figures 4, 8). We found that a more sophisticated training scheme and, for spiral waves, an increased NN capacity (5 layers of 64 neurons vs. 4 layers of 32 neurons, see Supplementary Table 2) were necessary for convergence in these large and complex problems. This more complex setup could, of course, have been used to solve the simpler 1D problems, through a trade-off between computational time and the convenience of a one-size-fits-all EP-PINNs approach.

It is in inverse mode, when estimating model properties from sparse measurements of V, that the EP-PINNs framework showcases its usefulness. Parameter estimation is an important topic in EP, as it is essential for both the personalisation of models and for the understanding of the effect of pathology and drugs on APs. Although parameter estimation in EP has been extensively discussed as a means of reducing the uncertainty associated with current biophysical models (34), NNs had not yet been used for dedicated model parameter estimation in EP. This contrasts with the more common use of NNs as efficient solvers of EP systems (in a similar fashion to the EP-PINNs forward mode in the current study) (15, 35, 36).

Using EP-PINNs, we were able to estimate, with different degrees of accuracy, three different biophysical parameters, each controlling, in an almost uncoupled manner, different observable properties of the system: APD (through b), excitability (through a) and conduction velocity (through D). In both 1D and 2D (homogeneous and heterogeneous) problems, the network was highly successful at estimating b, but struggled with a, especially when estimating it in tandem with D. This is likely to reflect a dependency between EP-PINNs' inverse mode effectiveness and the solution type that is probed experimentally. Indeed, when compared to b, the experimental inputs to the network (VGT(x→,t)), depend little on a providing the initial stimulus is supra-threshold. An exception could have been the spiral wave scenario, whose properties (e.g., the distance between spiral arms) depend strongly on model parameters such as a. Model parameters are more strongly coupled in the properties of spiral solutions, however, explaining the consistently lower accuracy of EP-PINNs' estimates in this scenario (Figures 7a,d). Parameter estimation in the spiral wave scenario may be improved when EP-PINNs are trained in longer time series, in which the spiral wave tip samples more of the spatial domain.

These issues underlie the difficulties of finding a single experimental design that allows for simultaneous accurate estimation of several EP parameters. A solution for this problem may be the training of PINNs using data from the same system acquired in different experimental conditions, perhaps by training separate NNs in parallel with a combined loss function, in an analogous manner to architecture B in this study (see Figure 2).

EP-PINNs were additionally able to identify heterogeneities in D across the 2D domain (see Figure 8), by estimating D in a separate parallel NN which shared terms of the loss function (L_data, see Equation 5) with the main NN. Characterising heterogeneities in EP parameters from electrical measurements is an interesting problem from a clinical point of view, as arrhythmias such as AF are often accompanied by heterogeneous changes in EP properties. Of these, regions of dense fibrosis [often modelled as areas with a reduced D (10)] are a promising candidate for personalised ablation sites (37), which may increase the overall efficacy of these procedures. EP-PINNs are thus well placed to help locate these putative ablation sites by identifying spatial heterogeneities in EP parameters such as D.

The inputs to the current EP-PINNs implementation, V(x→,t) are not, however, measurable clinically. Future work will modify EP-PINNs to instead perform parameter inference from extracellular electrical potentials, ϕ_e, which are regularly measured during clinical procedures using contact electrodes. ϕ_e can be interpreted as a weighted spatial convolution of the ∇→.(D∇→V) term in Equation (2) (38, 39), making the identification of localised EP changes more difficult. To take this into account, EP-PINNs designed for ϕ_e analysis may benefit from a move away from the current fully-connected architecture to incorporate, for example, convolutional layers.

An important point for future clinical applications of EP-PINNs is its ability to generalise beyond the details of the setup it is trained on. We showed that the proposed EP-PINNs implementation is model-agnostic, as it was able to perform robust parameter inference on in silico data generated by a much more complex atrial EP model (8) than the 2-variable ventricular one used in its loss function. In particular, EP-PINNs were able to correctly identify the decrease in APD (manifest as an increase in b) that is associated with increasing degrees of AF-induced remodelling in this model (see Figure 6).

In contrast to most previous PINN studies (17, 22, 28, 33, 40), we complemented the in silico studies with an assessment of the EP-PINNs performance on experimental biological data. Despite requiring a proportionally higher amount of training data than in silico experiments, EP-PINNs were able to cope well with the noise and artefacts unavoidably present in experimental data to identify the effect on APD of two different drugs: an I_Kr blocker and an I_CaL blocker. As for the data generated by a different mathematical model (Figure 6), EP-PINNs coped well with differences between the experimental data and the Aliev-Panfilov model, namely in resting membrane potential (see Figure 9). This ability to generalise well to data with different characteristics could be due to the use of PDE as a soft constraint (a term in the loss function) in the EP-PINNs framework, as well as the lack of assumptions about the distributions from which data come from.

As demonstrated for the in silico tests, the EP-PINNs framework can easily be extended to simultaneously infer the effect of drugs on more than one EP parameter, which may be useful for the characterisation and safety assessments of anti-arrhythmic drugs. These applications may further benefit from the training of EP-PINNs on more complex biophysical models, to obtain a more fine-grained characterisation of potential pharmacological (or pathological) effects.

The current study aims to demonstrate the potential of PINNs within EP, as an initial necessary step towards clinical applications of this method. As such, we only trained EP-PINNs using one comparatively simple EP model, in 1D and 2D scenarios and for the estimation of a small number of EP parameters. We additionally did not test EP-PINNs in the chaotic or pseudo-chaotic scenarios of spiral wave break-up (9), which may be relevant for some arrhythmias.

Generalisations of the proposed framework to 2D/3D geometries representative of cardiac chambers, anisotropic conditions and more detailed EP biophysical models can be achieved by further increasing the capacity of the deployed EP-PINNs, with a concurrent increase in computational resources. However, promising and less resource-intensive applications for EP-PINNs may be the characterisation of pharmacological effects on AP or the identification of heterogeneities in EP properties from EGM signals.