Reconstruction of excitation waves from mechanical deformation using physics-informed neural networks

Generation of synthetic dataForward modelThe synthetic data set used in this work was generated by modeling electrical activation waves which propagate through a 2D elastic medium. The physical equations for the excitation wave, coupling and mechanical deformation as well as the parameter values follow the work of a previous publication29. First, modified Aliev-Panfilov kinetics30 are used to model the excitation wave dynamics:$$\begin{aligned} \frac{\partial v}{\partial t} = \Delta v – kv(v-a)(v-1) – vw \end{aligned}$$
(1)
$$\begin{aligned} \frac{\partial w}{\partial t} = \epsilon (v)(kv-w) \end{aligned}$$
(2)
where v and w are the normalized transmembrane voltage and recovery variable, respectively. The step function \(\epsilon (v)\) is equal to 0.1 for \(v\ge a\) and 1.0 for \(v<a\) as in previous work31. Parameter values were set to \(k=8\), \(a=0.05\). This model has dimensionless space units (s.u.) and time units (t.u.). The action potential duration for \(90\,\%\) repolarization (\(APD_{90}\)) was measured in simulations to be \(6.8\,\)t.u.. Scaling this to a human ventricular \(APD_{90}\) of 270 ms32,33, 1 time unit in the model corresponds to 40 ms. Similarly, the measured conduction velocity for a one-dimensional isolated traveling pulse was 1.6 s.u./t.u. By comparing this to a tabulated value of wave propagation in human ventricle (0.68 mm/ms34), we conclude that one space unit in our model corresponds to 17 mm. In a next modeling step, the transmembrane voltage v induces active tension \(T_a\) in the medium according to the phenomenological ordinary differential equation35:$$\begin{aligned} \frac{\partial T_a}{\partial t} = \epsilon (v)(k_Tv-T_a), \end{aligned}$$
(3)
resulting in a similar \(T_a\) wave, slightly delayed in time with respect to the voltage. We used \(k_T=1.5\). Lastly, for the deformation resulting from the active tension \(T_a\), the linear and isotropic Navier–Cauchy equations for mechanical equilibrium in an elastic medium are used36$$\begin{aligned} (\lambda +\mu )\vec {\nabla }(\vec {\nabla }\cdot \vec {U}) + \mu \Delta \vec {U} = – \vec {\nabla } T_a \end{aligned}$$
(4)
Here, the passive elastic stresses, involving the material constants \(\lambda \) and \(\mu \) and spatial derivatives of the deformation vector \(\vec {U}\), are balanced by the active tension gradient. We used the parameter values \(\lambda =1.25\), \(\mu =1\). Note that since we work with linear elasticity and dimensionless active tension here, only the ratio \(\lambda / \mu \) is fixed by our parameter choice.For simplicity, we used a 2D square domain with no-displacement boundaries. This choice corresponds to the experimental setting where a left-ventricular tissue wedge is subtended in a rigid square frame. We are aware that our model is inherently limited, since we adopt linear elasticity, isotropy of excitation and contraction, and a simplified geometry and boundary conditions. Nonetheless, the set of Eqs. (1–4) provides a minimal model for the study of excitation-contraction coupling in the cardiac muscle.SimulationsAll simulations were carried out on a homogeneous and isotropic 2D square domain of size 100 s.u. for a time duration equal to 60 t.u.. Due to the absence of mechano-electrical feedback and the small displacement regime, the same static domain was applicable for every timestep and the physical equations could be solved sequentially. This introduces only a small error, but allows for great flexibility in simulation choices. The Aliev-Panfilov equations for the evolution of the transmembrane voltage (Eqs. 1, 2) were solved using the !ithildin package37, a finite difference time-stepper in C++ with MPI-parallellization, which was developed in our group and made publicly available. Simulations were executed on a 2D regular grid with a spatial resolution of 0.5 s.u.(8.5 mm), while integrating in time using forward Euler stepping at a temporal resolution of 0.05 t.u.(2 ms). Standard Neumann boundary conditions for the transmembrane voltage were applied to close the equations. We adopted 2 different initial conditions. One corresponded to 2 focal pulses at \(t=0\) and \(t=20\), initiated by inserting external voltage (\(v=1\)) in a circular region of radius \(r=1\) with centers located at (0, 0) and (75, 25), respectively. The second simulation involved a stable rotating spiral initiated via a S1-S2 protocol. Every 1 t.u. the transmembrane voltage was stored for all spatial grid points. Active tension was developed according to Eq. (3). The ODE was solved for every spatial point on a reduced grid (resolution 2.5 x 2.5 s.u.) with the use of python package scipy38 and its integration module integrate.odeint. The resulting active tension values were stored every 2 t.u.. Finally, we made use of the publicly available finite-element package Fenicsx39 to solve the Navier-Cauchy equations (Eq. 4). The regular grid points of the active tension solution were used as vertices for the quadrilateral elements of the finite element mesh, while the function space for both the active tension and deformation vector were set to the standard \(P_1\) linear Lagrange element. No-displacement boundary conditions were applied to all four sides. Final data consisted out of active tension and x- and y-components of the deformation vector on a 40 x 40 spatial grid (resolution 2.5 s.u.) for 30 time steps (resolution 2 t.u.).An example of a synthetic data set obtained by forward modeling is shown in Fig. 1. Two focal sources were initiated by stimulating a circular area of radius 1, differing in time and location. The active tension profile closely follows the voltage wave with only a small time delay and a slightly broader wave profile. As we are primarily focused in this study on recovering the accurate activation pattern, we will from here on out consider only the active tension \(T_a\) as the variable of interest and investigate how it can be reconstructed from its resulting deformation. Both deformation components show large values in regions around the \(T_a\)-wave, while also exhibiting non-negligible values further away, due to the global spatial coupling in the Navier–Cauchy equations (Eq. 4).Figure 1Time series of synthetic data generated as outlined in the Methods section. Two focal sources were initiated by inserting external voltage inside a circular region. The resulting transmembrane voltage wave produces a similar active tension wave, slightly delayed and broadened. The active tension then instantaneously induces deformation in the elastic medium.Reduction of data qualityWe tested the robustness of the developed method against the effect of noise, lower spatial resolution, and the absence of out-of-plane deformation components for tri-planar sampling. To this purpose, we processed the synthetic data in the following manner. We applied additive Gaussian noise to the simulated displacements \(U_{x,0}\) and \(U_{y,0}\). The noise level is taken proportional to the maximal absolute value of displacement in the x- and y-directions over the full spatiotemporal domain, i.e. \(U_{x,max} = {\textrm max}(U_{x,0})\), \(U_{y,max} = {\textrm max}(U_{y,0})\). Then, we have for every measurement point:$$\begin{aligned} U_{x} = U_{x,0} + \eta _{x}\ p\ U_{x,max} \end{aligned}$$
(5)
$$\begin{aligned} U_{y} = U_{y,0} + \eta _{y}\ p\ U_{y,max} \end{aligned}$$
(6)
Here, all \(\eta _{x}, \eta _{y}\) are sampled for every point in space and time from the standard normal distribution. The positive parameter p is the noise level of choice and was varied between \(0\%\) and \(50\%\) in steps of \(10\%\). Signal-to-noise ratios (SNR) were then calculated by dividing the variance of the original deformation values (signal) by the variance of the generated noise. This is done for both deformation components at every time step and afterwards averaged to obtain one SNR value for a specific data set. For a second group of data sets, we reduced the spatial resolution using a moving spatial average, since spatial resolution in an actual ultrasound recording depends on the scan sequence and hardware. Deformation components were averaged over \(n^2\) data points in a n x n square resulting in images used for inversion which have 40/n by 40/n pixels. n was set to 2, 4 and 8. Our third type of non-ideality was inspired by tri-planar acquisition of echocardiagraphic data in which three planes are simultaneously recorded in order to approximate the full 3D tissue displacement with high frame rate 2D methods40. The planes all share a common direction but are rotated around this axis (e.g. all planes are apical views where the common direction is the left ventricle longitudinal coordinate). To simulate this in our 2D data set, 40 points were sampled on three lines which were rotated by \({60}^{\circ }\) with respect to each other. This was done by sampling each line with 40 equally distanced points and for every point selecting the data coordinate and deformation values closest to it. An example of the data slices is shown in Fig. 8. As a true recording only contains information about the displacement within the viewing plane, the deformation vector \(\vec {U}\) was projected onto that line, resulting in one projected deformation scalar per point on the line, yielding only 3\(\times 40 = 120\) data points per time frame to be used as input data for the inversion.Quality assessment of the reconstructionTo evaluate the quality of reconstruction, we did not adapt a pixel-based metric such as the root-mean square error or the L2-difference between the ground truth and the reconstructed image. The underlying reason is that, in view of localizing the electrical sources of arrhythmia, we are primarily interested in the accurate shape of the activation wave rather than the absolute values. Additionally, there is only information present about the gradient of the active tension field through the Navier-Cauchy equations, such that the absolute baseline value of \(T_a\) cannot be determined from deformation alone. Therefore, we adopted in this work a well-established image quality measure called the Feature Similarity (fsim) index41. The fsim index is calculated from the phase congruence, a dimensionless quantity of the significance of a local structure, as well as the image gradient magnitude, see details in the original paper41. By this choice, we will have no undesired effects from higher background values and our results are evaluated with a focus on the global shape. Fsim values shown in this work for time series are calculated separately for every frame and then averaged over the time series to obtain a single score between 0 and 1 that covers the quality of reconstruction.Inversion methodsOur working hypothesis is that the inversion of sparse ultrasound measurements can be translated to a large optimization problem, where the unknown is the full spatiotemporal active tension field. Such large-scale optimization problem can be tackled by physics-informed neural networks, in which the neural networks that were previously designed for supervised machine learning tasks, are used as flexible non-linear function approximators. They are equipped with efficient optimization schemes and minimize loss functions that incorporate explicit physical knowledge, thus combining data-driven and model-driven methodologies.Navier–Cauchy physics-informed neural network (NC-PINN)To reconstruct waves of active tension from deformation data, a physics-informed neural network was implemented in TensorFlow42 (version 2.10.0), based on the point-wise formulation with a fully-connected, dense structure20 (\(N_l\) hidden layers of each \(N_n\) nodes), see Fig. 2. The neural network is used as a flexible function approximator, which maps spatial and temporal coordinates (x, y, t) to both deformation components (\(U_x, U_y \)) and active tension (\(T_a\)). We kept the \(\vec {U}\)- and \(T_a\)-prediction in the same network (in contrast to other work24), such that the identified spatiotemporal patterns in the deformation space can easily be shared with the active tension prediction. PINN optimization involves the minimization of a custom loss function that will ensure a match with the measured data, obey physical laws via differential equations and respect certain boundary conditions. The relative importance of these contributions is tuned via weighing factors \(\alpha _*\), to be discussed below. The following loss function \(\mathscr {L}\) was chosen in this work:$$\begin{aligned} \mathscr {L}= & {} \mathscr {L}_{data} + \mathscr {L}_{nc} + \mathscr {L}_{bound} \end{aligned}$$
(7)
$$\begin{aligned} \mathscr {L}_{data}= & {} \frac{\alpha _{data}}{N_d}\sum \limits _{i=1}^{N_d} || \vec {U}(x_i) – \hat{\vec {U}}_i||^2 \end{aligned}$$
(8)
$$\begin{aligned} \mathscr {L}_{nc}= & {} \frac{\alpha _{nc}}{N_c}\sum \limits _{i=1}^{N_c} || (\lambda +\mu )\vec {\nabla }(\vec {\nabla }\cdot \vec {U}(x_i)) + \mu \Delta \vec {U}(x_i) + \vec {\nabla } T_a(x_i)||^2 \end{aligned}$$
(9)
$$\begin{aligned} \mathscr {L}_{bound}= & {} \frac{\alpha _{bound} }{N_b}\sum \limits _{i=1}^{N_b} || \vec {U}(x_i) – \vec {0} || ^2 \end{aligned}$$
(10)
The first term \(\mathscr {L}_{data}\) calculates the mean squared error between the predicted deformation components \(\vec {U}(x_i)\) and the deformation data that are available \(\hat{\vec {U}}_i\). This calculation is done for \(N_d\) randomly sampled data points. Dividing by the number of points converts the sum to averages, such that the value is independent of the chosen number of sampled points. The second term \(\mathscr {L}_{nc}\) imposes that the output satisfies, at least locally, the Navier–Cauchy equations and is thus formulated in terms of spatial derivatives of the deformation and active tension as well as the mechanical parameters of the equation, \(\lambda \) and \(\mu \). Both parameters are assumed to be known in this study. Derivatives are conveniently calculated through automatic differentiation in the same manner as gradients for the NN parameters \(\theta \) are obtained. Because the physical equations are supposed to hold in the whole spatiotemporal domain, \(N_c\) collocation points are sampled using the Latin hypercube design. Lastly, the term \(\mathscr {L}_{bound}\) translates the no-displacement boundary conditions. Similarly, \(N_b\) collocation points are sampled through the Latin hypercube design, while setting here either the x- or y-coordinate equal to 0 or L with L the size of the spatial domain in order to obtain boundary points.Once the total loss function \(\mathscr {L}\) is calculated, network parameters \(\theta \) can be updated based on the gradients of the minimization of \(\mathscr {L}\) with respect to \(\theta \). This whole process is done for N iterations. We have set \(N=10^{4}\) and chose per iteration a subsample, often referred to as a mini-batch, of 64 points meaning \(N_d=N_c=N_b=64\). The optimization is done with the Adam optimizer of TensorFlow and its default parameters. A schematic overview of one iteration is shown in Fig. 2. To facilitate the optimization, spatial and temporal coordinates were first normalized to the range \(\left[ 0,1\right] \), while the active tension output was put through the TensorFlow sigmoid function and multiplied by a maximum value before calculating the loss function. This approach is possible as \(T_a\) will always be non-negative and an upper limit can be well estimated, similar to the magnitude of velocity in other studies24. Lastly, the hyperbolic tangent function was used as the differentiable activation function and Glorot initialization from a uniform distribution has been applied to all NN parameters before the start of optimization.Figure 2Schematic overview of one iteration of the optimization process in the NC-PINN. The neural network consist of a densely connected architecture of \(N_l\) hidden layers with \(N_n\) nodes each. Inputs are spatiotemporal coordinates (x, y, t), outputs are deformation components \(U_x\), \(U_y\) and active tension \(T_a\). Colours indicate the points (coordinates) used for each loss term. The network is updated after the total loss is minimized with respect to the network parameters \(\theta \) (weights and biases of the network).Additional regularization using wave propagation constraint (EIK-PINN)In the case of poor data quality, we propose to use additional physical information. For simple wave patterns without repeated activity, the isotropic and homogeneous eikonal equation43 can be used. This viewpoint relates the local activation time \(\tau (\vec {r})\) to the local conduction velocity \(v(\vec {r})\):$$\begin{aligned} \Vert \vec {\nabla }\tau (\vec {r}) \Vert = \frac{1}{v(\vec {r})}. \end{aligned}$$
(11)
Although detailed measurements may allow to estimate a local conduction velocity24, \(v(\vec {r})\) was taken to be constant in space in this work. Moreover, curvature effects on the propagation speed44,45 are neglected. A schematic overview of the method combining the elastic and wave propagation regularization is provided in Fig. 3b. First, the deformation data was fed into the NC-PINN as outlined in Fig. 2, to reconstruct a first estimate of the active tension, denoted \({\hat{T}}_a\). For every spatial coordinate on a fine regular grid, the initial local activation time \({\hat{\tau }}\) was then found by searching for the time when \({\hat{T}}_a(t)\) attains its maximum over the full time domain. Only points for which \(T_a({\hat{\tau }}) \ge T_{a,max}/3\), with \(T_{a,max}\) a predefined constant, were kept. Finally, these \({\hat{\tau }}\) values were used as input to a second PINN, with wave propagation regularization, see Fig. 3a. Since this PINN serves as an approximator to the function \(\tau (\vec {r})\), it only takes x and y as input, not time. The loss function for this network was:$$\begin{aligned} \mathscr {L}’= & {} \mathscr {L}_{data-\tau } + \mathscr {L}_{eik} \end{aligned}$$
(12)
$$\begin{aligned} \mathscr {L}_{data-\tau }= & {} \frac{\alpha _{data-\tau }}{N_d}\sum \limits _{i}^{N_d} (\tau (x_i) – {\hat{\tau }}_i)^2 \end{aligned}$$
(13)
$$\begin{aligned} \mathscr {L}_{eik}= & {} \frac{\alpha _{eik}}{N_c}\sum \limits _{i}^{N_c} \left( \Vert \vec {\nabla }\tau (x_i) \Vert -\frac{1}{v}\right) ^2. \end{aligned}$$
(14)
The data term consists of the pre-calculated \({\hat{\tau }}\) values, while the eikonal term supplements additional information in order to inter- or extrapolate activation times or to overwrite potential errors from the previous NC-solution. To return to the \(T_a(x,y,t)\) space, we make use of a predefined \(T_a(t)\) shape to approximate the time evolution of typical active tension progression. In this work we have assumed that the tension developed over time follows a Gaussian curve, a simplifying but sufficient function, of given width \(\sigma \):$$\begin{aligned} T_a(t;\tau ) = T_{a,max} \exp \left( \frac{(t-\tau )^2}{2\sigma ^2} \right) \end{aligned}$$
(15)
with constant parameters \(T_{a,max}\) and \(\sigma \). We have chosen \(v=1.5\), \(\sigma \) = 2.5 and \(T_{a,max} = 0.42\), as approximate values of the measurements from the data set (\(v_{exp}=1.6\), \(\sigma _{exp} = 2.3\), \(T_{a,max,exp} = 0.42\)). In applications, these values will have to be either estimated based on patient-specific data, on literature values, or in the case of the conduction velocity, could be predicted simultaneously with the activation times. The EIK-PINN could in principle be coupled directly to the NC-PINN by only estimating the \(\tau \)-field and calculating \(T_a\) based on Eq. (15), eliminating the need for the \(\tau \)-data term. Preliminary results with this simultaneous approach, however, showed highly unstable training processes leading to erroneous predictions. Therefore, we designed and report here the two-step approach of Fig. 3b.Figure 3(a) Overview of the EIK-PINN to enhance the reconstruction of the local activation time (LAT) \(\tau \). Colours indicate the points (coordinates) used for each loss term. (b) Schematic overview of the two-step optimization (NC+EIK) in case both elastic and wave propagation constraints are used. The standard NC-PINN is used to convert deformation \(\vec {U}\) into a first estimate of active tension \({\hat{T}}_a\). From this field, the LAT is calculated, and optimized subject to the eikonal relation (Eq. 11) via a second PINN, depicted in panel (a). The emerging \(\tau \) can either inter- or extrapolate the solution or overwrite initially difficult regions, after which it is converted back to \(T_a\) via Eq. (15).