Robust data assimilation with noise: Applications to cardiac dynamics

Weather forecasting has improved dramatically over the last few decades due not only to improvements in understanding and modeling atmospheric physics as well as increases in computing power but also due to the development of specialized data assimilation algorithms incorporating real-world measurements to form a more complete picture of meso-scale weather patterns.⁴ A similarly difficult forecasting problem is found in cardiac arrhythmias, in which spatiotemporally complex patterns interact with both anatomical features and self-organized dynamical structures such as phase singularities and scroll wave filaments.^5–7 These dynamical regimes are distinguished by sustained electrical turbulence across and within cardiac tissues, in which small changes in the configuration of organizing features may lead to dramatically different outcomes over short times.⁸ This difficulty for accurate forecasting is compounded by sensitivity to details of the present configuration of electrical signals in tissue and the unobservability of the details of these complex states. For example, in the context of ex vivo experiments on animal hearts using optical mapping measurements, only a minimal depth below the surface of the tissue is available to be measured,⁹ and then typically only the transmembrane potential, leaving the distribution of the remaining physiologically relevant ionic concentrations (e.g., calcium, sodium, etc.) and other state variables undetermined.

The cardiac forecaster follows in the steps of the weather forecaster, aiming to reconstruct the present state of the system from often sparse and noisy observations. Such observations are irreducibly sparse and noisy in certain contexts—notably, in clinical applications.¹⁰ Such techniques are often referred to as “state reconstruction” under the banner of data assimilation, which takes uncertain observations of a state and approximates the underlying dynamics that lead to those observations; techniques that exploit the synchronization of loosely coupled systems generally [i.e., our Local Ensemble Transform Kalman Filter (LETKF) ensemble states and truth model state] exist under the same name with similar goals of state and model parameter estimation.^11,12 Various approaches to data assimilation appear in the literature (c.f. Ref. 13 for a review of available methodologies) but we shall focus on those which utilize an ensemble of simulated states to forecast, broadly Kalman filter (KF) type methods and specifically the Local Ensemble Transform Kalman Filter (LETKF).¹⁴ Our use of LETKF is in part motivated by the structure of the model system while weather system models typically include difficult linear terms with modest nonlinearities, cardiac systems typically couple modest linear terms with stiff and sometimes discontinuous nonlinearities, making the linearization of the nonlinear models especially difficult. Data assimilation systems that utilize adjoint tangent evolution (e.g., “weak” 4D-Var¹⁵) introduce significant overhead through the recording of the forward model for the evaluation of the adjoint evolution. Likewise, we utilize a sequential filter approach to data assimilation (as opposed to “smoother” methods) with consideration for real-time applications in which the future configuration of the system is strictly inaccessible. In these methods, the predicted state of the system is represented by the mean of a model-generated ensemble which, when observed, approximates the observations of the true state with some uncertainty.

A particular issue with ensemble data assimilation methods is that the variance of an ensemble, or ensemble spread, can collapse while the ensemble mean disagrees with the features of the truth–convergence to an inaccurate reconstruction. As the ensemble spread is a proxy for the uncertainty of the state reconstruction, this signals greater certainty in an inaccurate state estimate. Such ensemble collapse is more likely when the underlying ensemble model is strongly nonlinear and the leading instabilities of the system are under-represented in a finite ensemble of much smaller dimension than the original state space. The ensemble spread is thus a subject of scrutiny in these methods: increasing the variation within an ensemble may improve the accuracy of the reconstruction by effectively sowing doubt and allowing corrections arising from new observations to more strongly affect the ensemble members, and, indirectly, the ensemble mean. The maintenance of ensemble spread is handled by inflationary processes either within the assimilation process or within the evaluation of the model. Within the assimilation process, the ensemble variance may be maintained by a scaling of the deviations from the ensemble mean—multiplicative inflation—and by changing the configuration of the deviations from the ensemble mean by the addition of new, typically stochastic, terms—additive inflation. The maintenance of ensemble spread through the model equations themselves, and the effect on the ensemble mean error as a result, is the primary focus of this investigation.

Forecasting generally requires knowing not only the present state of the system but also a model of the dynamical relations underpinning the evolution of the system. Cardiac researchers, thankfully, have spent years developing both generic and specialized mathematical models of various cardiac tissues.¹⁶ However, even in principle, models are unable to capture the exact relations within these complicated physiological systems, and thus model error is endemic and should be addressed in attempts to forecast future states. Model error is by no means a new concern; weather forecasting has repeatedly investigated methods for correcting the bias of model systems,^17–19 with bias-removal approaches using statistics generated from a precomputed history in both state-dependent and state-independent correction schemes having found particular success.^20,21 Approaches that model the propagation of model error explicitly^22–24 have also shown improvements relative to otherwise uncorrected assimilation algorithms. Ensemble-based data assimilation algorithms generally address model error using the existing features of the process (the ensembles) through inflation schemes,^25–27 while still others have combined bias removal with inflation schemes in ensemble methods to great effect.²⁸ In this investigation, we consider the inclusion of stochasticity in data assimilation processes strictly through model definitions to account for systemic model error and to reduce state reconstruction error for sparsely observed dynamics.

The use of stochastic models of deterministic processes in constructing probabilistic forecasts goes back at least 60 years.^29,30 Stochastic parameterizations have been included that perturb deterministic models when the effect of specific model terms is ambiguous,³¹ while further revisions incorporating additional physical motivations and stochastically parameterized effects have found improvements, especially in weather forecasting.³² Still others have considered explicit multi-model approaches^33–35 or schemes with parametric variation.³⁶ Some have argued explicitly that systemic model errors are due in part to the coupling of large-scale deterministic truncation of the flow and stochastic modeling of sub-grid processes.³⁷ Comparisons of multi-model and multi-parameter techniques in the service of seasonal ensemble forecasts performed similarly to other, state-of-the-art, approaches.³⁸ Stochastic parameter dynamics, an “artificial noise evolution model,” have also been shown to improve assimilation results with a neuronal model using mismatched parameters when included in the assimilated state.³⁹ 

In this work, we focus on traditionally stochastic model descriptions and on stochastic parameterization methods for cardiac systems, with an eye toward managing ensemble spread through model modifications.

This work is laid out as follows. In Sec. II, we present a mathematical description of the cardiac excitation model used in this work and a foundation of the particular data-assimilation method we utilize in the forecasting of this cardiac model. We further describe the methods utilized for numerical experiments, including the mathematical model of cardiac excitation and stochastic methods, in Sec. III. Section IV shows results from experiments in small and large cardiac model systems with different observation procedures. Section V presents some final discussion and future directions.

Forecasting the future state of a system requires both a model of the dynamics and a method of estimating the current state of the system from observations; the two interact through the generation of a state estimate and the correction of the state estimate according to observations of the true state. In the context of ensemble data assimilation methods, the reconstructed states are described by ensemble means—an average over several realizations of the dynamical state from multiple model estimates. We discuss details of both the model and the assimilation framework below, before finally discussing stochasticity with regards to the model construction.

The Fenton–Karma model describes the dynamics of the transmembrane voltage $u$⁠, sodium-mediated gating $v$⁠, and calcium-mediated gating $w$ at the level of homogeneous tissue.⁴⁰ The model equations are expressed using a partial differential equation (PDE) for the transmembrane potential $u(t,x)$⁠, with $x$ denoting a position in one or three dimensions, coupled to ordinary differential equations (ODEs) governing the dynamics of gating variables $v(t,x)$ and $w(t,x)$ mediating the calcium and sodium currents. The model equations are given by

The PDE is accompanied by boundary conditions on $u(t,x)$⁠, which in this case are either “no-flux” (⁠$n^⋅∇xu(t,x)=0$⁠) or periodic. The Heaviside function is denoted by $Θ(x)$⁠, such that $Θ(x)=0$ for $x<0$⁠, $Θ(x)=1/2$ for $x=0$⁠, and $Θ(x)=1$ for $x>0$⁠. The local relaxation time of $v$ depends on the local value of $u$⁠: $τv−(u)=Θ(u−uv)τv1−+Θ(uv−u)τv2−$⁠.

The ionic currents correspond to a fast inward current gated by $v$ representing sodium, a slow inward current gated by $w$ representing calcium, and a slow outward current representing potassium that is mediated by self-interaction with the transmembrane potential $u$⁠, $Jion(u,v,w)=Jso(u)+Jfi(u,v)+Jsi(u,w)$⁠; these currents are described by the following nonlinear functions, each contributing to the morphology of excitations in space and time, or “action potentials,”

The dynamics of the gating variables are semi-linear, with short transitions between growing and decaying dynamics, and with spatial coordination controlled by the interaction with the transmembrane potential, $u(t,x)$⁠, which is continuous in space (⁠$C2$⁠) and time (⁠$C1$⁠).

For a fixed spatial resolution, we define the configuration of the state combining the transmembrane potential with the gating variables, $x=[u(x),v(x),w(x)]$⁠, where $x$ regularly and discretely indexes over the spatial domain. We define the model $M$ as an operator that maps $x(t)→x(t+Δt)=Mx(t)$⁠, where $M$ solves the initial-boundary value problem defined by (1) with appropriate boundary conditions and parameter definitions over the prescribed time interval $Δt$⁠.

Data assimilation (DA) is a discipline of related mathematical and computational techniques that aim to systematically incorporate observations into a mathematical model of the system. This procedure requires quantification of the uncertainty in both the background state and the observations. A number of algorithmic approaches to DA exist and are used primarily by the weather forecasting community, including variational approaches^41,42 and those with history derived from the Kalman filter (KF).⁴³ While the former category benefits from explicit modeling of the linearized dynamics, it propagates observation information both forward and backward in time, making online applications more difficult. The Kalman filter explicitly models the state and covariance update from the previous state and covariance estimate, and the lineage of DA algorithms that derive from the KF have been applied to especially large models. Broadly speaking, the Kalman filter involves an analysis and background state estimate, where the analysis step is a maximum likelihood estimate relating the background estimate of the current state to observations of the state estimate. One such essential extension of the KF to larger systems is the ensemble Kalman filter (EnKF),^44–46 which treats the derivation of the error covariance relation from the sampling of an ensemble of state vectors itself rather than solving the expensive covariance update equations. The ensemble transform Kalman filter (ETKF)⁴⁷ is a further specialization of the KF algorithm that combines the ensemble approach of EnKF with an error-subspace transformation, which more efficiently updates the ensemble state vectors and covariances.

The ETKF algorithm, in the context of twin experiments relevant to the present study, begins with a “truth” sequence $xo(t)∈Rm$ that describes the state of a model system over time. The application of a (potentially nonlinear) noisy observation operator generates the truth observations, $yo(t)=H(xo(t))+ϵo(t)$⁠, where $ϵo(t)∼N(0,σo2)$ and $yo(t)∈Rl$ with, $l≪m$⁠. A “background” ensemble of states generated by the dynamics of a model system, $xb=[x1b,x2b,…,xkb]$⁠, is used to generate an “analysis” ensemble of states $xa=[x1a,x2a,…,xka]$⁠, where the analysis ensemble mean $x¯a$ estimates the current state. The analysis ensemble mean $x¯a$ minimizes the Kalman filter functional,

where $Pb=(k−1)−1∑i=1k(xib−x¯b)(xib−x¯b)⊤$ is the background covariance, $R$ is the observation operator covariance matrix, and $x¯b=k−1∑i=1kxib$ is the mean of the background ensemble. The $m×m$ background covariance $Pb$ is never formed explicitly; instead, the transformation from background ensemble to analysis ensemble is formed in the $k≪m$-dimensional ensemble space. To make such updates explicit, we relate the mean $y¯b$ and mean deviation $Yb$ of the observations of the background ensemble members, $yib=H(xib)$⁠, providing a basis with which to relate the ensemble and observation spaces, $Yb=yb−y¯b$⁠. In turn, we define the analysis ensemble covariance in the background ensemble basis, $P~a=[(k−1)I+(Yb)⊤R−1Yb]−1$⁠, where $P~a∈Rk×k$⁠. The weighting of each mean deviation of the background ensemble determines the transformation from $xb$ to $xa$⁠; for clarity, we consider the update of the mean (⁠$x¯b→x¯a$⁠) and of the mean deviations (⁠$Xb→Xa$⁠) separately. The analysis ensemble mean $x¯a$ is determined by the mean weightings in the background ensemble basis, $w¯a=P~a(Yb)⊤R−1(yo−y¯b)$⁠, such that

What remains is to update the individual mean deviations of the analysis ensemble and form $xa$⁠. The weighting of the mean deviations is likewise expressed in the ensemble space, where we define the weights using the symmetric square root, $Wa=[(k−1)P~a]1/2$⁠. From the recombination of the analysis weights $wa=w¯a+Wa$⁠, it follows that

where we have defined the abstract (nonlinear) ETKF operator $K1$⁠.

The algorithmic innovation relevant to the study and assimilation of large systems is the incorporation of covariance and observation localization, designating the local ensemble transform Kalman filter (LETKF).¹⁴ Localization of the observation influence and covariance of the ensemble exploits the locally low-dimensional dynamics of real physical systems and enforces the decorrelation of long-range dynamics. The LETKF makes the assimilation of large dynamical models computationally tractable, splitting, and parallelizing the $m$-dimensional system into a set of $l$-dimensional matrix systems, where $l≪m$ denotes the dimension of the local subspace. We refer the reader to Ref. 14 for a complete and detailed expression of the ensemble update equations deriving from (2), especially the localization transformations in concert with (3) and (4); for brevity, we assert that the application of ETKF to each element of a state vector with appropriately localized observations captures the relevant concept.

The method proceeds to correct the mismatch between the true state and the analysis ensemble mean so long as the Kalman update is capable of correcting when presented with new observations. When the ensemble consensus is inaccurate and the ensemble cannot be updated to account for new observations, the method may converge to an inaccurate state reconstruction due to systematic underestimation of the uncertainty of the state estimate. An approach to preventing this breakdown is to maintain adequate ensemble variance through the perturbation of the state vectors relative to the ensemble mean. The ensemble variance may be maintained through the incorporation of multiplicative and additive inflation processes. Multiplicative inflation may be applied to the analysis covariance through multiplication by the desired inflation factor $ρ≥1$⁠, where $ρ=1$ corresponds to no inflation. Multiplicative inflation, therefore, is incorporated entirely within the generation of the analysis ensemble from the background ensemble; for this reason we include multiplicative inflation within the Kalman update operation, $Kρ$⁠.

Just as the LETKF is only able to correct errors that lie in the span of the background ensemble deviations from the mean, multiplicative inflation is typically effective for ensembles large enough to contain the leading instabilities of the model and the deviation from the truth state. Additive inflation, on the other hand, perturbs the analysis ensemble used in the generation of the updated background ensemble, although the precise perturbation method varies.^48,49 In this way, additive inflation includes new perturbation directions in the initial conditions of the background ensemble.

Finally, we combine the model, inflation, and assimilation steps into the update of the background and analysis ensembles from time $t=nΔt$ to time $t=(n+1)Δt$⁠,

where the deterministic evolution $M$⁠, additive inflation $Aα$⁠, and the Kalman filter assimilation process with multiplicative inflation $Kρ$⁠, all act sequentially.

The form of (5) is suggestive: we may consider instead a stochastic evolution operator $Mσ$ that incorporates continuous stochastic forcing of the ensemble member dynamics rather than the dynamics-less additive inflation given by $Aα$⁠. That is, beginning from (1), we append a stochastic forcing term to the dynamical model that defines a SDE,

where the terms inside $[…]$ represent the deterministic model (1), and the stochastic forcing term, $g(u)$⁠, is to be determined. The noise vector $dWt$ denotes a vector of independent and identically distributed (i.i.d.) Gaussian noise processes evaluated over all space and time with fixed variance and mean. In practice, this is the generation of a random number from a prescribed normal distribution $N(0,σu2)$ for each numerical degree of freedom, which is modified by the action of $g(u)$⁠. This amounts to reformulating (5) as

interpreting the additive inflation as a stochastic process within the model, rather than an additional process post-Kalman filter. So, where the deterministic evolution operator $M=M0$ solves (1), the stochastic evolution operator $Mσ$ solves (6). This interpretation presents an opportunity to motivate the inflation—to drive the ensemble members according to the underlying instabilities of the dynamical system—and to tune the functional form of the inflation in a physically motivated way through the construction of specialized $g(u)$⁠.

We shall investigate the effects of some different functional forms for the stochastic term of (6) in the assimilation process. We initially consider forcing with a diagonal stochastic term, which drives all state variables using noise processes with spatially uniform i.i.d. parameters,

which perturbs the evolution of every element of each state vector with an independent noise process. In the limit of short assimilation intervals (with the assimilation evaluated after every time step $Δt$⁠), this term is equivalent to the application of $Aα$⁠, with $α=σuΔt$⁠.

We proceed from the forcing defined by (8) noise by considering the time scales associated with the model. First we split the stochastic dynamics along the same fast–slow splitting as the deterministic model, forming two new stochastic terms that consist of switching the stochastic forcing to the voltage channel or the gating channels. The voltage channel and gating channels forcing terms are written as

respectively.

The approach outlined above is motivated solely by the fast–slow dynamics of the deterministic system. By testing the efficacy of (8), (9), and (10) noise, we can determine effective channels for the reduction of ensemble mean error during the assimilation process. Further tweaking the form of $g(u)$ may improve the efficacy of these methods by targeting the largest error or the smallest ensemble spread. In addition to reinforcing the time-scale splitting of the underlying dynamics, the forcing terms (9) and (10) perturb the observed and unobserved state variables, respectively. Therefore, which components of the state vector are observed should be carefully considered. In this study, we have chosen to observe the transmembrane potential variable $u(t,x)$⁠, as the row of $∂xMx$ corresponding to perturbations in $u$ maintains full rank in a single time step; other choices (e.g., $v$ and $w$⁠) require multiple time steps for state derivatives to achieve full rank. It is expected that, while these derivatives are not used in the assimilation method, any approximation of the linear neighborhood of the state through other means should similarly reflect the sparsity of the true derivative, and thus may suffer if information about unobserved states is propagated only over longer time intervals and through the dynamics of a coupled variable, e.g., observing only $v$ may lead to poor convergence due to information about the state of $w$ being inaccessible except through its effect on $u$⁠, which interacts with $v$ directly. Observations of the transmembrane potential also reflect the observations most readily accessible in experiments.

An alternative to the SDE formulation of the dynamics that has a direct and simple analogy to classical inflation methods is the stochastic selection of model parameters (SMPs) for each ensemble member during the construction of the background ensemble, $xb$⁠. This method may be written analogously to (5) and (7) as

where $Mi,n′$ is a model parameterization of (1) selected stochastically by perturbing the model parameters, indexed by the ensemble index $i=1…k$ and the assimilation step $n$⁠. That is, we define $σp$ as a vector of the same dimension of $po$—the true model parameter Vector—and for each ensemble member model invocation we draw a randomized parameter vector $p$ such that

When $σpi=0$⁠, then $pi=p0i$⁠, and when $σpi≠0$⁠, then $pi$ is drawn from a Gaussian distribution centered on $p0i$ with variance $σpi$⁠. The parameters are selected once for each ensemble member on each assimilation interval—i.e., each time the model is invoked, each parameter is selected from its respective distribution.

The SMP method is analogous to a restricted application of additive inflation. It acts once (not continuously) per ensemble member model evaluation, and discretely and discontinuously impacts the ensemble member dynamics. Whereas additive inflation updates the ensemble members $xi→(xi)′$ prior to time-evolution by model $M$⁠, this technique updates the model $M′$ while preserving the continuity of the input state vector $xi$⁠. In this way, where additive inflation and the SDE method act to introduce new directions to the ensemble corresponding to relevant instabilities of the underlying dynamics (whether physically targeted or not), the SMP method attempts to model uncertainties in the parameterization of the model directly.

In principle, the parameterization of the model may be only partially stochastically determined, e.g., so that only the set of parametric time scales or parameters related to the excitation threshold are perturbed in the generation of $M′$⁠. Likewise, for a model with underdetermined parameters, one may select an appropriate variance to capture the true dynamics with mis-specified ensemble models. We will consider two SMP rules; the first selects all the explicit time scales of the deterministic model (the $τ∗$⁠) stochastically,

called the SMP-$τ$ method. The second selects the parameters associated with the excitation threshold,

called the SMP-$c$ method. The disparity of scales across the parameters of the model presents an issue to the scaling of the noise; therefore, both methods use the deterministic values of the parameter to scale the per-parameter variance appropriately; otherwise, it is likely that, e.g., for small $τ∗¯$ the SMP-$τ$ method would cause a divergence in the model simulation by selecting $τ∗≈0$⁠. The SMP-$τ$ method impacts all the time scales of the dynamics, similarly to additive inflation and (8) noise, and so we might expect similar performance in the accuracy of the reconstruction. The excitation parameter selection (SMP-$c$⁠) addresses only the fastest dynamics of the system, specifically those related to the wavefront, which may have similar short-term effects to SMP-$τ$⁠, but may exhibit longer term drift in the error of the state reconstruction.

We choose these two SMP rules analogously to the SDE terms; we consider a temporally broad stochasticity affecting all state variables and time scales of the system in SMP-$τ$ [analogously to (8) noise], and specialize to the fastest time scales of the system and excitation thresholds using SMP-$c$ [analogously to a further targeted (9) noise]. Parametric stochasticity that affects the slowest state variables or dynamics are no further specialized than (10) noise and so we omit such a stochastic parameter rule.

The efficacy of an assimilation scheme in the absence of a truth sequence, as in an experimental setting when only the observations are available, is a subtle question. “Twin experiments,” in which the truth state is accessible, make the assessment of an assimilation scheme much more straightforward. We define several quantities used to assess the assimilation reconstruction.

The assimilation reconstruction is compared to the truth sequence $xo(t)$ by computing the difference from the analysis or background ensemble mean $x¯a,b(t)$⁠,

such that $eia,b(t)$ describes the error in the analysis or background ensemble averages at state vector index $i$ and time $t$⁠. Due to the online nature of the reconstruction process, it is useful to consider specifically the root-mean-squared (RMS) error as a singular measure of the reconstruction error over a desirable subset of the state described by the index mask $q$⁠,

such that $|q|$ is the length of the indices spanned by $q$⁠. Additional measures that incorporate information about the system by summing along only a subset of physical dimensions (e.g., $q={z}$⁠) or only over a subset of state variables (e.g., $q={v,w}$⁠), or some combination (e.g., $q={z,u}$⁠), are also useful measures in specific contexts. It is sometimes useful to consider the time-averaged $eRMSa,b(t)$ for sufficiently long experiments,

where the assimilation interval $ti−ti−1$ is shorter than the time scale of asymptotic dynamics of the system, while $tN$ is much longer than both the assimilation interval and the longest time scales of the underlying dynamics.

Additionally, we define the ensemble spread,

which is the pointwise variance of the ensemble state vectors at time $t$⁠, where $i=1,…,k$ indexes the ensemble vectors. Similarly to $eRMSa,b(t)$⁠, we may compute the RMS of the ensemble spread as with (15), privileging certain state variables or a subset of physical dimensions, to form an appropriate $sprdRMS(t)$⁠,

This quantity can be treated as an approximation of the variance of the ensemble and a proxy for the error of the state reconstruction when the truth sequence is inaccessible, i.e., during the assimilation of experimental observations. Ideally, the ensemble spread closely predicts the ensemble mean error of the truth state reconstruction.

Throughout this study, we shall consider the ensemble mean error and ensemble spread over all the spatial dimensions of the $u$-variable, i.e., $q={u,x,…}$⁠, unless otherwise specified. These metrics thus privilege the most readily observed variable in experimental settings.

In this section, we describe our numerical models and assimilation experiment programs. The model systems include a one-dimensional system in a periodic dynamical regime and a three-dimensional system with fiber rotation generating anisotropic diffusion along and across the local fiber orientation exhibiting sustained spiral turbulence.⁵⁰ The LETKF 1D implementation is coded in Matlab, while the 3D implementation is customized based on the FORTRAN code written by Takemasa Miyoshi (https://code.google.com/p/miyoshi/) and run using Comet at the San Diego Supercomputing Center through XSEDE.⁵¹ We specify for each case the details of the assimilation process, including the ensemble initialization schemes, inflation techniques, and localization parameters.

The one-dimensional model is evaluated on a domain discretized with $N=560$ grid points and grid size $δx=0.025cm$⁠. The spatial coupling term is evaluated with a sparse second-order centered-difference stencil corresponding to periodic boundary conditions and spatial homogeneity, i.e., $∇xD(x)⋅∇xu(x)≡0$⁠. The model is time-stepped using a forward-Euler update with fixed time step $δt=0.05ms$⁠.

The truth model develops alternans, which is a periodic variation of the width⁵² and amplitude⁵³ of an excitation wave as it propagates around the domain, eventually converging to a stable relative periodic orbit solution, c.f. Fig. 1. This relaxation to the stable alternans dynamics⁸ is the truth state sequence, $xo(t)$⁠, against which the efficacy of the assimilation reconstruction is compared.

We introduce model error into the system by assigning a mismatch in the diffusion constant of the truth and ensemble models. This mismatch is intended to investigate some of the errors associated with the calculation and use of the diffusion coefficient in experimental settings. Diffusive properties of cardiac tissue can vary not only intramurally but also across the surface of the heart from the apex to the base.⁵⁴ For the generation of the truth sequence, we assign the diffusion constant to be $D||′=(0.9)−2D||$⁠, where $D||$ is the diffusion constant of the ensemble model(s), otherwise leaving all parameters the same (c.f. Table I). Likewise, the parameter sets listed in Table I are chosen to reproduce particular dynamical regimes in the corresponding geometries: alternans, which is fundamental to dynamical cardiac modeling, in 1D and scroll wave turbulence, relevant for ventricular fibrillation, in 3D. The specific parameter sets in 1D and 3D experiments were adapted from existing parameterizations of the Fenton–Karma model.⁵⁵ The difference from the ensemble model diffusion to the truth model diffusion is a comparatively large $19%$⁠, equivalent to a domain size change of $10%$⁠, from $14.0cm$ (ensemble) to $15.56cm$ (truth). The spatial resolution and the diffusion constant are unaffected by any of the inflation methods utilized and are not included in the assimilation system nor in the observations of the ensemble members. In this way, all forecast ensembles that originate from exactly reproduced truth states are guaranteed to diverge under uncorrected model evolution.

We initialize the background ensemble using a random sampling of the history of the truth state, i.e., $xib(t=0)=xo(ti)+ξi$⁠, where $ξi∼N(0,σo2)$ and $ti$ is selected from the recent history, $ti∼U[−40,0)ms$⁠. Numerical experiments modifying the ensemble initialization indicate that this method generates both larger initial ensemble mean error and larger ensemble spread compared to, e.g., selecting the $k$ most recent samples from the truth history for an ensemble of size $k$⁠, on average, and thus presents a more difficult assimilation task than some alternative initialization methods. We use an assimilation window of $5ms$ and run each experiment up to $t=2000ms$ with an ensemble size of $k=6$⁠.

Smaller ensembles are desirable for several reasons; foremost is the speed of execution and reduction of computational resources. However, small ensembles typically suffer in ensemble data assimilation schemes due to the need to properly model the relevant range of instabilities in a physical system; small ensembles systematically underestimate the spread and provide an overconfident prediction. In the interest of faster iteration on methodologies, it appears that more effective use of a small set of ensemble members is more desirable than simply increasing ensemble sizes. One such way is to use “fair scores” of assimilation methods,⁵⁶ such as the continuous ranked probability score (CRPS),⁵⁷ which may be trivially transformed to account for expected improvements due to larger ensembles for methodology comparisons.⁵⁸ In any case, small ensembles present a challenge for many ensemble-based data assimilation techniques, and making more effective use of each ensemble member presents an opportunity.

We observe the state at $35$ sites (every $0.875cm$ or $16$ grid points) across the domain, uniformly on the state variable $u$⁠. The observation localization parameter is $hloc=2$⁠, which defines the decaying influence of each observation on the covariance. The covariance localization for each observation is truncated at $⌊210/3hloc⌋=7$ grid points or $0.175cm$⁠, so that uncorrected physical regions of the domain are small (two grid points or $0.05cm$⁠) and uniformly spread. Each observation includes added noise generated from $N(0.0,0.052)$⁠, which distorts the voltage measurement [c.f. Figs. 2(a)–2(d)], especially in quiescent regions. The observation noise is modeled with zero bias because, for all but unreliably small observation dimension, LETKF manages observation bias automatically by considering the difference between the observations of the truth state and the observations of the ensemble member states. This level of noise [the half-height-width of $N(0.0,0.052)$] is representative of experimental observations using optical mapping (⁠$8%--10%/100mV$⁠).⁵⁹ 

We begin by assessing the utility of classical inflation methods—multiplicative and additive inflation—controlled by the inflation parameters $ρ$ and $α$⁠. We treat these assimilation experiments as baselines against which the stochastic inflation schemes are compared. In these experiments, additive inflation adds perturbations to the ensemble from a randomly selected subset of sequential differences of the truth state generated during “spin-up”; each perturbation vector is scaled by the additive inflation parameter $α$⁠. We consider multiplicative and additive inflation individually before considering their combination. These experiments correspond to nine assimilation experiments, with $ρ−1,α∈{0.00,0.05,0.10}$⁠, compared over time and to the observation error.

With the baseline established, we consider the effects of the SDE terms on the assimilation. We consider the form of the SDE terms, comparing (8) to (10) and (9) noise terms for a fixed variance $σu$⁠. Likewise, we consider the effect of changing SDE variance on the reconstruction error for a demonstrably effective SDE term, namely, (8) noise.

Then, we investigate the efficacy of the SMP methods, SMP-$τ$ and SMP-$c$⁠. In particular, we compare the SMP-$τ$ and SMP-$c$ methods in the absence of any other inflation processes with fixed variances. We then consider the interaction between SMP-$τ$ and multiplicative and additive inflation schemes, which is roughly the limit of prescribed sampling for the optimization of the ensemble mean error without more sophisticated algorithms.

Finally, we treat the selection of inflation scheme parameters $(ρ,α,σu,σp)$ as a stochastic optimization problem. Rather than a single monolithic problem, we consider four in sequence: the optimization of $e¯RMS$ over $(ρ,α)$⁠; over $(ρ,σu)$ for each noise term; over $(ρ,σp)$ for SMP-$τ$⁠; and finally $(σu,σp)$ for each noise term with SMP-$τ$⁠. These results are the first truly nonlinear investigation of the inflation schemes in this study as the limitation of “small” inflation parameters is relaxed.

The three-dimensional model is run on a numerical grid of size $200×200×50$ points, with grid spacing $δx=0.02cm$⁠, a domain size of $4×4×1cm3$ (similar to the thickness of a slab of left ventricular muscle in humans⁶⁰). The diffusion of the transmembrane potential is computed using a centered-difference approximation with anisotropy due to fiber orientation, with $D∥=0.001$ and $D⊥=0.0002cm2/ms$ along and across the local fiber direction, respectively. The fiber orientation varies with depth,⁶¹ changing by an angle of $π/3rad$ (⁠$60°$⁠) over the $1cm$ range of the $z$ coordinate, where the rotation of the fiber orientation is explicitly $θ(z)=(π/3)(z−1/2)$⁠, so that the diffusion is spatially heterogeneous. No-flux boundary conditions are applied to all domain boundaries and evaluated using the ghost-point method in concert with the finite-difference stencil. The model is time-stepped using an operator-splitting approach with $δt=0.025ms$⁠; the deterministic portion of the gating variable dynamics (⁠$v$⁠, $w$⁠) is updated with a first-order exponential time-differencing method.⁶² This is done strictly for convenience for this model, as the semilinear dynamics of $v$ and $w$ encourage an analytical stepping scheme not limited by the stability limits of an explicit scheme. The nonlinear and diffusive updates for the transmembrane potential $u$ are evaluated using a simple forward-Euler scheme. The stochastic terms are handled separately for all state variables, also incorporated using an Euler–Maruyama method.

The observations are localized to the $z=0$ and $z=1$ levels (grid indices $1$ and $50$⁠), taken every $3$ grid points in the $x−y$ plane, corresponding to a spacing of $0.06cm$⁠. The localization scale is $hloc=6$ grid points, giving a covariance localization truncation radius just shy of $22$ grid points (with an isotropic Gaussian taper) or $0.438cm$⁠, such that the center the domain is only updated indirectly—by the model evolution—during the assimilation process. For this reason, we expect the assimilation efficacy to vary with the $z$-depth of the state, and likewise for dynamics in which information flow is primarily outward from the interior of the domain toward the $z$-boundaries to generate intermittent spikes in the overall error of the state reconstruction.

The base cycle period in this system is 125–130 ms, so the configuration of the system dramatically changes between each time slice in Fig. 3. The 3D model system sustains scroll wave turbulence over time; the excitation waves are themselves organized about the scroll wave filaments, which connect two counter-rotating waves at $t=0$ s and $z=1$ cm, to two counter-rotating waves at $t=2$ s and $z=0$ cm, cf. Fig. 3(b). These organizing filament dynamics happens primarily in the interior of the domain. For the observation organization for this system—restricted to the $z=0,1$ cm planes—the scroll wave filament configuration presents a significant challenge, as the bulk of the organizing features of the state appear far from the observed regions and the most significant changes propagate from the inner part of the domain to the boundaries.

The ensemble is composed of $k=20$ members generated from the history of the truth model, selected sequentially from the 20 assimilation times prior to the initial truth state, i.e., $xi(t=0)=xo(t−iΔt)$⁠, for $i=1,…,k$⁠. No additional noise is added, distinguishing the initialization of the 3D ensemble from that of the 1D ensemble. As this ensemble procedure spans a longer history, exactly reproduces configurations accessible to the truth model, and includes closely correlated configurations of the state in the initial ensemble, this initialization procedure yields an ensemble, which more likely contains the dominant features of the true state at the initial time than the method that is used in the 1D experiments, which only contains the original truth state on average and in the limit of large ensembles.

We consider the baseline experiment in 3D, with small multiplicative inflation (⁠$ρ=1.01$⁠) and no additive inflation (⁠$α=0$⁠). We then compare the inclusion of noise in the evolution model, the inclusion of stochastically perturbed timescales, and then the combination of SDE and SMP-$τ$ methods on the ensemble mean error and spread. In each experiment, we present the integrated ensemble mean error and ensemble spread over both the $x$ and $y$ coordinates and over either the $t$ or $z$ coordinates.

We apply stochastic assimilation techniques to the Fenton–Karma cardiac model in both one and three spatial dimensions. The one-dimensional application emphasizes the utility of the assimilation for addressing sources of endemic model error that cannot be reduced by typical model parameter modifications. First, we consider assimilation using only classical multiplicative and additive inflation methods. We then compare the classical inflation methods to stochastic differential equation forms of ensemble models and models with stochastically perturbed system parameters. These include SDE forms corresponding to (8), (9), and (10) noise, as well as the SMP-$τ$ and SMP-$c$ methods.

The three-dimensional application addresses the utility of the stochastic methods with sparse observations of relatively large systems. As in the low-dimensional experiments, we compare classical inflation methods (multiplicative inflation only) to the application of stochastic models [SDEs corresponding to (8) only] and SMP-$τ$⁠. In concert, the one-dimensional and three-dimensional experiments explore some defining issues related to data assimilation of cardiac signals.

We assimilate the 1D model dynamics under LETKF first with no inflation processes (⁠$ρ=1$⁠, $α=0$⁠) and then a combination of multiplicative and additive inflation methods (⁠$ρ={1.05,1.10}$⁠, $α={0.05,0.10}$⁠). Figures 4(a) and 4(b) depict the background and analysis ensemble mean RMS errors compared to the observation error for the state reconstruction. Without inflation (⁠$ρ=1$⁠, $α=0$⁠; blue), the LETKF is not able to accurately reconstruct the truth state; the RMS error in both the background and analysis ensemble means decreases to $0.2$ before increasing to $O(1)$—comparable to the norm of the state variable $u(t,x)$⁠, representing a wholly inaccurate ensemble mean. The ensemble spreads shown in Figs. 4(c) and 4(d) are likewise indicative of a failure to reconstruct the state and decay to $O(10−3)$ after $400ms$⁠. As expected, multiplicative inflation $ρ>1$ (orange) is effective at maintaining a larger ensemble spread over time [c.f. Figs. 4(c) and 4(d)], although the ensemble mean error is not improved dramatically. The incorporation of additive inflation (⁠$α>0$⁠; green) improves both the maintenance of ensemble spread and the RMS error of the ensemble mean; however, the RMS error is still larger than the observation error for nearly all times. The inclusion of both multiplicative and additive inflation (red) makes dramatic improvements to both the RMS ensemble errors and the ensemble spreads, although only a marginal improvement over additive inflation alone. These runs serve as baselines against which we will compare the improvements to the ensemble mean error and ensemble spread using stochastic inflation methods.

Figure 5 depicts the result of assimilation in the absence of traditional multiplicative and additive inflation methods (⁠$ρ=1$⁠, $α=0$⁠) but with strong stochastic dynamics in the ensemble model (⁠$σu=0.1$⁠, $σp=0$⁠) using the (8), (9), and (10) noise terms. It is here that we begin to see the dynamical fast–slow splitting of the model interact with the assimilation; stochasticity is a “fast” process—the noise is accumulated on the scale of the time step, which by construction is shorter than any of the dynamical time scales of the system. The fast stochastic forcing methods [(8) (blue) and (9) (green) noise terms] improve on the ensemble mean error compared to the best multiplicative and additive inflation offers [c.f. Figs. 4(a), 4(b), and 5(a)], although the ensemble spread is comparable [c.f. Figs. 4(c), 4(d), and Fig. 5(b)]. The accumulation of stochastic noise in slow channels, i.e., (10) noise term (orange), is competitive with multiplicative inflation alone [compare Figs. 4(a) and 4(b) with Fig. 5(a)]. Alternatively, we can see periodic decreases in the ensemble mean error 5(a) and periodic increases in the ensemble spread 5(b) associated with the (10) noise experiment. These behaviors correspond to the collapse of the action potential width due to alternans and thus the momentary dominance of the gating variable dynamics; the effect is subtle, apparent in the spread when the error is small and apparent in the error when the spread is small. Likewise, while the error of the (10) noise result exhibits less apparent variability than the errors resulting from the other noise terms; this is due primarily to the logarithmic scaling of the error and spread, cf. Figs. 5(a) and 5(b), and secondarily to the measurement of the $u$-variable error while the stochastic term primarily affects the gating variables.

A consequence of the stochastic model approach to the reduction of ensemble mean error is the convergence of the background and analysis ensembles. Equation (7) suggests an explanation: in the absence of additive and multiplicative inflation (⁠$A0=I$⁠, and $ρ=1$⁠), any explicit ensemble inflation is completely localized to the evaluation of the model, $Mσ$⁠, and so any distinct improvement in ensemble spread or ensemble mean error (from background to analysis ensembles, $xb$ to $xa$⁠) must be due solely to the action of the Kalman filter, $K1$⁠. In this scenario, the Kalman filter $K1$ is contracting from $xb$ to $xa$ (as we should expect when all the ensemble members are perturbed to fit the same set of observations). Likewise, the stochastic forcing in the ensemble model is the sole generator of ensemble divergence; when the assimilation interval increases, the stochastic forcing terms generate more ensemble spread between assimilation steps. For short assimilation intervals (such as those effective for Kalman filter based assimilation approaches), we might expect the coalescence of the background and analysis ensembles, while for longer assimilation intervals, we should expect larger differences between the background and analysis ensembles.

Figure 6 compares the SMP-$τ$ and SMP-$c$ stochastic inflation methods for small values of $σp$⁠. These experiments indicate that SMP-$τ$ is much more effective than SMP-$c$ for reducing the ensemble mean error, Fig. 6(a), at the same variance parameter $σp$⁠. This is not entirely surprising: while SMP-$c$ is directly associated with the excitation front parameters and the error of the reconstruction, it affects only the transmembrane potential directly and is thus a significantly smaller state perturbation than SMP-$τ$⁠. In another sense, the specialization of SMP-$c$ leaves a larger space of perturbations to the ensemble inaccessible than SMP-$τ$⁠, which affects all state variables through the modification of the model time scales. Indeed, SMP-$τ$ with $σp=0.05$ (blue) yields a smaller error than SMP-$c$ with $σp=0.05$ (green) or $σp=0.10$ (red). The spread due to the SMP-$τ$ and SMP-$c$ methods [Fig. 6(b)] is generally lower than that generated by any of the fast SDE methods; this may be due to the difference in the continuous forcing of the fast SDE methods and the forcing on the time scale of the assimilation interval due to the SMP-$τ$ and SMP-$c$ methods.

For sufficiently small model and ensemble sizes, we may treat the selection of inflation parameters as a stochastic optimization problem. We define the assimilation error measure to be minimized as the root mean square of the analysis ensemble error over time, (16), and pass the assimilation parameters $ρ$⁠, $α$⁠, $σu$⁠, and $σp$⁠, to the assimilation process for a combination SMP-$τ$ and SDE method using either (8), (9), or (10). We use the stochastic minimization function patternsearch available in the Matlab optimization toolbox—this method requires no derivatives of the objective function and instead samples the function sequentially, following an adaptively refined gradient descent from successfully sampled starting points.

We optimize the value of $eRMSa$ as defined by (16) for the analysis ensemble $xa$⁠, by sampling one or more of $ρ$⁠, $α$⁠, $σu$⁠, and $σp$⁠. The optimization cannot give a globally optimal set of parameters for this model and method, but it can inform about the relative importance of the inflationary parameter strengths and their utility for a given dynamical regime, e.g., alternans. The form of the objective measure is likewise motivated not by the absolute minimum that $eRMSa(t)$ achieves at some time, but by sustained accurate reconstruction of the truth state. Alternatively, we may privilege the background ensemble mean error $e¯RMSb$ for purely short-term predictive matters (although evidence presented thus far indicates that the results would be largely unaffected), or compute some alternative scoring measure, e.g., continuous ranked probability score (CRPS),⁵⁷ which attempts to measure the confidence of the analysis and background ensembles without recourse to the truth sequence, i.e., $CRPS(xb,x¯a)$⁠.⁵⁶ 

We proceed by considering (i) an optimization of only the multiplicative and additive inflation parameters (⁠$ρ$⁠, $α$⁠); (ii) an optimization of the combination of multiplicative and SDE parameters (⁠$ρ$⁠, $σu$⁠); (iii) an optimization of the multiplicative and SMP-$τ$ parameters (⁠$ρ$⁠, $σp$⁠); and (iv) an optimization of the stochastic parameters without any multiplicative or additive inflation (⁠$σu$⁠, $σp$⁠). In each case where $σu>0$⁠, we likewise optimize using each of the stochastic forcing terms: (8), (9), and (10). Each optimization procedure is initialized from the no-inflation case $(ρ=1,α=0,σu=0,σp=0)$⁠, with the appropriate subset of parameters permitted to update according to the optimization.

Table II describes the parameters computed from these optimization problems. These results largely echo those presented in Figs. 4–6; given the source of model error, strong parameter noise using SMP-$τ$ is effective at reducing the ensemble mean error regardless of additional multiplicative inflation or stochastic state dynamics. Likewise, stochastic state dynamics targeting fast timescales (8) and (9) is effective in concert with some small amount of multiplicative inflation, although (10) is unable to reduce the ensemble mean error, in particular, that associated with the deviation in the wavefront propagation speed due to mismatched diffusion. The comparison of stochastic methods (ii)–(iv) with combination multiplicative and additive inflation (i) is strongly in favor of the former—the latter is unable to attain sufficiently low ensemble mean error during the optimization process.

Figure 7 shows the analysis ensemble mean errors for the optimized parameters listed in Table II, both (a) over time and (b) as a time-average (for the given context of model error driven by parametric mismatch, and initial condition uncertainties). Several features should be noted; case (i), the classical inflation techniques in tandem (blue) are unable to reduce the analysis ensemble mean error to below the level of the observation error, $σobs=0.05$⁠, even for large amounts of multiplicative and additive inflation. Additionally, case (ii) is infrequently successful, with only the (8) noise (orange) outperforming (on average) the classical methods—notably while using the smallest multiplicative inflation factor and smallest noise variance of the set— while (9) and (10) noise (green and red, respectively) do not. The analysis ensemble mean error for the optimized case (iv) is, however, uniformly successful (8), (9), and (10 noise; brown, pink, and gray, respectively). Each experiment successfully reduces the time-average error to below the observation error, and each optimization uses similar values for $σu$ and $σp$⁠. The remaining experiment, corresponding to the optimization of $(ρ,σp)$ for case (iii) achieves the lowest ensemble mean error (purple), though it requires several more assimilation steps to reduce below the observation error level, so the distribution of ensemble mean errors for case (iii) is mostly comparable to those of case (iv).

These results suggest a preferred sequence of parametric tuning for the reduction of the ensemble mean error. Let us take the case of (8) noise in Fig. 7; beginning from the default $(ρ=1,α=0,σu=0,σp=0)$⁠, we find that multiplicative and additive inflation (blue) reduces the reconstruction error compared to the default, but further increases are found through a combination of multiplicative inflation and SDE noise (orange). Further reductions are made using SDE noise in combination with stochastic parameterization of the ensemble model(s) (brown), reducing the reconstruction error below the observation error threshold. At each step, this parameter tuning improves the time-averaged ensemble mean error. The story is only superficially different for the (9) and (10) noise terms. Finally, we may further reduce the time-average ensemble mean error by considering multiplicative inflation in tandem with SMP-$τ$ (purple). However, in addition to increasing the time-to-observation-error-threshold, this final reduction in reconstruction error is distinguishable when the truth sequence is already known, as the combination of SDE model noise and SMP-$τ$ already suffices to bring the reconstruction below the observation error threshold.

Run-to-run variation of the efficacy of the stochastic inflation processes is of concern when the distribution of the ensemble mean errors exhibit “long tails,” potentially due to resonant effects between the nonlinearity of the model dynamics and the stochastic forcing. The efficacy of the optimization process suggests that the temporal variance of the ensemble mean RMS error is representative of the run-to-run variation. We have additionally accumulated each optimization case to generate distributions of the ensemble mean RMS error and ensemble RMS spread (see Fig. 11 in the  Appendix) for the (8) noise term and SMP-$τ$⁠, when appropriate. These distributions systematically over-sample parameter values near the default case compared to strict run-to-run variation sampling, but they still depict good convergence toward the optimal cases. A consideration of the variance computed about the maximum probability of the ensemble mean RMS error gives the variation of the stochastic realizations about the optimal parameter values.

Due to the computational cost in running long assimilation schemes in three spatial dimensions, we have opted not to test variations in the stochastic parameters exhaustively. Instead, we demonstrate the effectiveness of the stochastic methods by choosing parameter values that are small (to avoid over-driving the more complex dynamics of the 3D model) and that were found to be effective in 1D, and by comparing these results to the deterministic LETKF method. Additionally, unlike the 1D simulations, in the 3D model, the diffusion is unchanged between the truth and the ensemble forecast model; in this way, we have only observational and initial condition uncertainties within a “perfect model” scenario.

Figure 8 depicts the ensemble mean error and ensemble spread of the background and analysis ensembles over time $t$⁠. The error is a reduction of the full error $e(t,x,y,z)$ by (15) to a function of time $eRMS(t)$ over the voltage variable channel; likewise, the spread is a reduction of $sprd(t,x,y,z)$ by (18) to $sprdRMS(t)$⁠. The baseline experiment with no stochastic forcing (⁠$σu=0$⁠, $σp=0$⁠; blue) shows that with only a small amount of multiplicative inflation (⁠$ρ=1.01$⁠) the method is periodically effective at limiting the growth of ensemble mean error with respect to both the uncertainty in the initial configuration of the state and uncertainties due to noisy observations. The ensemble mean errors, both background and analysis, are approximately equal and slowly grow over the course of the experiment in an oscillatory fashion. Likewise, the spread for both the background and analysis ensembles are nearly indistinguishable for the length of the experiment, with decay from the initial ensemble spread to roughly the observation error, $σo=0.05$⁠, for the baseline run. We note that the classical inflation methods are unable to reduce the ensemble mean error to levels comparable to the observation error. Likewise, while the ensemble spread approximately matches the observation error, it systematically underestimates the ensemble mean error.

The introduction of stochastic state dynamics using (8) noise (orange) has minimal effect on the distribution of analysis ensemble mean errors in time. Likewise, the ensemble spread is largely unaffected by the inclusion of a stochastic noise term, at least for the relatively small variance seen here, $σu=0.02$⁠. One explanation is that the ensemble already reflects the instabilities of the true solution dynamics, so that using noise to stimulate already existing modes in the ensemble is not an effective method of reducing the ensemble mean error, nor of increasing the spread of the ensemble. That is, the three-dimensional experiments suggest that SDE dynamics are effective for the maintenance of ensemble spread when the ensemble size is a limiting factor; stochasticity makes better effective use of smaller ensemble sizes. Alternatively, the more complex dynamics of the 3D simulation may be less sensitive to small perturbations afforded by the noise term, and more strongly constrained by the topological features of the state, i.e., the scroll wave filament.

Stochastic parameterization (green) is effective at reducing the background and analysis ensemble mean errors [c.f. 8(a) and 8(b)]. Likewise, the spreads of both ensembles are increased with the inclusion of SMP-$τ$⁠, better maintaining the initial spreads of both the background and analysis ensembles over time [c.f. 8(c) and 8(d)]. Using SMP-$τ$ in concert with (8) noise forcing (red) improves both the background and analysis ensemble mean error over time while maintaining the ensemble spread. However, while the reduction of ensemble mean error from SMP-$τ$ by the inclusion of (8) noise forcing is subtle, any improvement in the ensemble spread is negligible.

We may conclude, at a glance, that while stochastic inflation schemes are not as effective for reconstructing the turbulent scroll wave dynamics as they were for the alternans regime, they reduce the ensemble mean error and, importantly, they lead to more realistic estimates of the ensemble mean error using the ensemble spread.

The distribution of the background and analysis ensemble mean errors through the depth, $z$⁠, are shown in Figs. 9(a) and 9(b), respectively. These correspond to an RMS reduction over $(x,y)$ according to (15) and (18), followed by the temporal average over the run. It is in this spatial representation that the sensitivity to covariance localization is realized—the error is asymmetrically distributed through the bulk of the tissue with both $z=0,1$ corresponding to local minima. Our conclusions for the effects stochastic methods have on the morphology of the ensemble mean error distribution reflect those found for Figs. 8(a) and 8(b): SMP-$τ$ (green) achieves an occasionally dramatic reduction in ensemble mean errors, while the combination of SMP-$τ$ and (8) noise (red) yields further improvement. Figures 9(a) and 9(b) further reveal something that the time series obscures: (8) noise in isolation (orange) is ineffective at reducing the ensemble mean error and may transiently increase it enough to impact the time-average. Finally, we note that, on average, the maximum background ensemble mean error in $z$ for the combined stochastic method (red) is smaller than the minimum background ensemble mean error in $z$ for the baseline experiment (blue). That is, confidence about the region of tissue adjacent to the observation planes in the deterministic assimilation runs is now available in the completely unobserved interior of the domain.

Figure 10 depicts the distribution of the reconstruction error, $ua(t,x,y,z)−u(t,x,y,z)$⁠, sampled over $t$ and $z$⁠, where $ua(t,x,y,z)$ is the analysis ensemble mean from the LETKF with (a) no stochastic inflation and (b) $σu=0.02$ and $σp=0.02$ with (8) noise and SMP-$τ$⁠, compared with the truth field, $u(t,x,y,z)$⁠. These parameters correspond to those used for the generation of the (a) blue and (b) red line in Figs. 8 and 9. The spatial distribution of the error shows the improvement in the reconstruction error near the boundaries of $z$ that the assimilation manages. Within $0.5s$⁠, the $z=0,1cm$ planes reduce the error distribution both relative to the initial ensemble mean at $t=0$ and relative to the interior of the domain.

As the ensemble mean does not satisfy the model equations, some features of the dynamics (especially, in this instance, the scroll wave filaments) are unlikely to be reliably determined from the ensemble mean and are not guaranteed to correspond to the corresponding features of the truth state. Efforts to incorporate the emergent features of coherent structures into data assimilation approaches have shown promise in systems in which the predominant contribution to the model error arises from continuous symmetries, namely (special) Euclidean symmetry.⁶³ Likewise, when linear ensemble means are used to reconstruct nonlinear states, the reconstruction is potentially unphysical for such systems. Heuristic approaches based on a limiter for the smoothness of the reconstruction may help us to reduce the error of the reconstruction by optimizing the weights of the ensemble sum for smooth solutions.⁶⁴ Physically motivated constraints have also been included in weather-forecasting assimilation settings to similar effect.^65,66 Robust recovery of the filament structure may benefit from a synthesis of these programs.

In this study, we have investigated formulations of ensemble inflation techniques that extend the underlying numerical model using stochastic effects applied to models of cardiac excitation in one and three spatial dimensions. These stochastic effects take the form of an extension of the deterministic numerical model to an SDE—an inclusion of an explicit stochastic forcing term that interacts with features of the desired dynamical system—or choosing a perturbation in the parameter space of the model definition. We find that stochastic effects that accumulate on time scales comparable to the shortest in the underlying deterministic system are effective at reducing ensemble mean error and maintaining appropriate ensemble spread. Likewise, in 3D applications with sparse observations, especially using both stochastic model dynamics and a stochastically perturbed model, stochasticity may significantly reduce the time-averaged ensemble mean error and better align the ensemble mean error with the ensemble spread.

We have characterized SDE model effects based on their relationship to the underlying nonlinear ensemble model, in this case a phenomenological cardiac model, and the inherent slow-fast splitting of the deterministic system. We have shown that, in the absence of multiplicative or additive inflation capabilities in the LETKF assimilation method, stochastic modifications of the ensemble model may qualitatively reproduce their effects. More precisely, modifications that affect the model on shorter dynamical time scales can closely match the effects of additive inflation (⁠$α>0$⁠), while modifications that affect the longer dynamical time scales are effective at reproducing the effects observed using multiplicative inflation (⁠$ρ>1$⁠). Indeed, while this qualitative correspondence helps clarify the underlying mechanisms at play, we found that stochastic methods were better able to reduce the ensemble mean error and maintain ensemble spread than either multiplicative or additive inflation using standard techniques. Even when the stochastic method at hand was unable to accurately reconstruct the truth state [e.g., using the SDE description with (10) noise], it tended to fare better than similarly equipped multiplicative inflation techniques. Adding multiplicative inflation to the assimilation process in addition to the stochastic forcing yielded marginal effects in the ensemble mean error and the ensemble spread.

This work has focused exclusively on the effect stochastic ensemble inflation techniques may have on the accuracy of state reconstruction; any number of other relevant parameters in LETKF are also likely candidates for improving the efficacy of these methods and may (or may not) interact in subtle and problematic ways. In particular, we have found through experimentation that the inclusion of stochastic effects may increase the sensitivity of the reconstruction to the localization processes of LETKF. Specifically, increasing the coverage of the observation localization so that the entire domain is modified in the construction of the analysis ensemble leads to strong coupling between the stochasticity of the model and the Kalman filter, in turn increasing the efficacy of other inflation parameters, notably multiplicative inflation $ρ$⁠. Further, as a model-specific modification, these stochastic inflation methods primarily affect the background ensemble; improvements in the analysis ensemble are due to additive and multiplicative inflation in concert with the Kalman filter. However, because the stochastic inflation techniques are model-specific, adapting these methods to work in tandem with traditional ensemble inflation techniques is programmatically straightforward.

Some such model-agnostic techniques are relaxation to prior posterior (RTPP)⁶⁷ and relaxation to prior spread (RTPS)⁴⁹ schemes, which manage the spread of ensembles in LETKF and similar data assimilation methods. Initial tests using the stochastic modeling approaches presented in this work suggest that RTPP and RTPS may over-drive the stochastic effects of the system and transiently increase the reconstruction error. Work is ongoing to characterize the effective regions of these techniques and their impact on stochastic inflation processes such as the SMP-$τ$ and (8) noise forcing used in this study.

While this study has considered model error through a parametric mismatch of the diffusion constant (or perfect model parameterization in the 3D experiments), the conclusions are not limited to only diffusion mismatches. Previous studies^2,3 included experiments considering classical inflation methods for model error due to kinetic parameter mismatches between the truth and ensemble models, namely, $τd$⁠, and for the total fiber rotation in 3D. We have considered similar kinetic model errors (see the  Appendix) and found that our main conclusions continue to hold. However, preliminary results suggest that some forms of model error may be somewhat easier to account for in LETKF with only multiplicative and additive inflation than the diffusive mismatch alone and are thus far from a worst-case scenario. In this light, work to provide a theoretical basis for optimally managing model error scenarios in a simplified model is ongoing.

A continuation in this approach will consider stochastic parameterizations of more flexible models across discrete sets of pre-tuned parameters capable of reproducing more distinct dynamics. Such efforts attempt to “bridge the gap” between multi-model, multi-parameter, and stochastic model parameter ensembles. Additionally, given careful consideration of relevant subsets of model parameters, the assimilation process can be extended to include model parameters explicitly,⁶⁸ leading to tuning of these models to observed data simultaneously with the state reconstruction.

Initialization of the ensemble has the potential to affect the accuracy or convergence rate of the ensemble mean. We have shown previously² that the effects of the random selection of states for ensemble initialization tend to be short-lived without altering the time frame for convergence (see Figs. 1 and 2 in this reference). Our observations in the current scenarios are consistent with these results.

This work considers 1D and 3D cardiac model applications with consideration for ventricular applications. Atrial cardiac modeling frequently utilizes a 2D domain geometry, which may be considered as a thin limit of the 3D domain as the depth decreases to only a few cell layers. With sufficiently dense observation distributions in the $x−y$ plane as the 3D slab, we expect these methods should be effective for atrial assimilation applications, as the central difficulty of the 3D application (unobserved dynamics in the interior of the domain) is eliminated. However, significant sparsity or non-uniformity of observations, which sometimes could arise in clinical studies, could pose greater challenges.

We have used stochastic forcing to reduce the ensemble mean error and increase the ensemble spread of assimilated observations of deterministic models. However, experimental observations are occasionally noisy not only due to an uncertain observation process but also due to genuinely stochastic underlying truth dynamics. We have begun work considering the potential for stochastic inflation processes (SMP-$τ$ and SDE models) for assimilating stochastic dynamics, where the stochastic forcing of the ensemble and truth models may agree or disagree, and comparing their efficacy to classical inflation schemes.

M.J.H. and E.M.C. were supported by the NSF under Grant Nos. CMMI-2011280 and CMMI-1762803. F.H.F. was supported by the NSF under Grant No. CMMI-1762553. C.D.M., F.H.F., and E.M.C. were supported by the NIH under Grant No. 1R01HL143450-01. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-1548562. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) Comet at San Diego Supercomputing Center at UC San Diego through allocation TG-IBN050000N.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

In any methodology utilizing stochastic effects, a consideration of the variance of the objects of inquiry with respect to the specific realization of randomness is important. In the course of the optimization schemes presented in Fig. 7, we generate a significant history of stochastic realizations, especially realizations with small (⁠$O(10−4)$⁠) deviations from the determined optimal parameters. From this History, we generate distributions for the analysis ensemble mean error and ensemble spread expected of these methods.

The distributions of analysis ensemble mean errors and ensemble spreads in Fig. 11 offer a succinct interpretation: the stochastic methods are systematically preferred for higher-accuracy reconstruction of the state sequence over the conventional scheme. Similarly, while stochastic methods incorporating multiplicative inflation may manage lower expectation values for the ensemble mean RMS error, these methods tend to have wider variance than the fully stochastic scheme, which may lead to spontaneously large errors in a single realization. Likewise, while the variation near the most likely realization of the ensemble mean error is larger for the (semi-) stochastic methods, these methods simultaneously exhibit faster decay away from the highest-probability realization while the conventional method exhibits (relatively) long tails.

We have considered model error via diffusive mismatch assimilation in 1D and model error free assimilation in 3D. An exhaustive theory of model parametric mismatch and the interaction with stochastic inflation effects is beyond the scope of this manuscript, but we can provide some intuition from numerical experiments focusing on kinetic parameter mismatches. Toward this end, we have considered two additional scenarios for model error due to parameter mismatch between the truth and ensemble models in 1D.

The first scenario modifies $τd→0.27$ and $uc→0.10$ in addition to the diffusive mismatch (c.f. Table I). Increasing $τd$ lowers the excitability of the model and decreases the conduction velocity of the excitation wave, decreasing $uc$ lowers the excitation threshold (increasing excitability), and increasing $D||$ increases the conduction velocity of the wave. Under these parameter modifications, the truth model exhibits very similar dynamics to the ensemble model with very different model parameters.

Conventional inflation schemes are similarly effective at smaller values of $(ρ,α)$ for the first scenario compared to those in the main text. This suggests that the present combination of parameter mismatch is somewhat more readily assimilated than the diffusive mismatch alone, as we might expect from more similar state dynamics. However, when stochastic methods are used instead of conventional inflation schemes, the stochastic methods outperform the conventional multiplicative and additive inflation schemes at both decreasing the ensemble mean error and ensuring the ensemble spread better represents the true error (not shown). When used in concert with classical inflation methods, the stochastic inflation methods reduce the error of the reconstruction and improve the spread of the ensemble compared to the conventional methods alone, often by nearly an order of magnitude (c.f. Fig. 12). This scenario of parameter mismatch highlights the identifiability problem of sufficiently complex models for data assimilation approaches.

The second scenario follows previous model kinetic mismatches,² where the ensemble model and truth model have the same parameter values save for $τd→0.29$⁠. In this scenario, the truth model relaxes to a periodic wave solution with minimal alternans while the ensemble model generates significant alternans in line with Fig. 1. This scenario presents a significantly more difficult assimilation task, which is reflected in the analysis ensemble mean RMS error (c.f. Fig. 13). The SDE method with 8 noise reduces the ensemble mean error when used instead of the conventional inflation scheme (and the SMP-$τ$ method, which is rendered ineffective in this scenario) but is unable to sustain errors comparable to the observation error (not shown). Similarly, the SDE method with (8) noise reduces the ensemble mean error when used in concert with conventional inflation schemes, while maintaining the ensemble spread at appropriate levels. Efforts to develop an exhaustive theory of parametric mismatch sensitivities in the context of state reconstruction for cardiac models are ongoing.