CoMetaPNS: Continually Meta-learning Personalized Neural Surrogates for Cardiac Electrophysiology Simulations

Graphical Abstract CoMetaPNS: Continually Meta-learning Personalized Neural Surrogates for Cardiac Electrophysiology Simulations

arXiv:2606.07488v1 [cs.LG] 5 Jun 2026

Ryan Missel, Xiajun Jiang, Linwei Wang

Highlights CoMetaPNS: Continually Meta-learning Personalized Neural Surrogates for Cardiac Electrophysiology Simulations Ryan Missel, Xiajun Jiang, Linwei Wang • Personalized heart models are hard to scale due to tuning time and compute cost • Meta-learning quickly personalizes neural surrogates from low subjectspecific data • Clinical data shifts over time, but existing meta-models assume fixed training sets • We extend meta-learning to continual learning using reservoir and clustering tools • Our method adapts to new subjects without retraining and avoids forgetting past data

CoMetaPNS: Continually Meta-learning Personalized Neural Surrogates for Cardiac Electrophysiology Simulations Ryan Missela , Xiajun Jiangb , Linwei Wanga a

Golisano College of Computing and Information Sciences, Rochester Institute of Technology, Rochester, NY, USA b Department of Computer Science, Rowan University, Glassboro, NJ, USA c The University of Utah„ Salt Lake City, UT, USA

Abstract Personalized virtual heart simulations face challenges in model personalization and computational cost. While neural surrogates offer state-of-the-art solutions, they typically address either efficient personalization or training generalizable models. Recent work reframes this by learning the process of personalizing a surrogate using limited subject-specific context data, through few-shot generative modeling with set-conditioned surrogates and meta-learned amortized inference. These methods, however, assume a static and diverse training distribution with known task identifiers. When new data becomes available, they require costly retraining with all prior data to avoid catastrophic forgetting - a phenomena where the model forgets earlier tasks when trained on new ones. This is a major limitation in clinical settings where often unlabeled data arrives sequentially and full retraining is infeasible. This paper presents a new continual meta-learning framework to achieve personalized neural surrogates able to not only continually integrate information but also identify whether incoming data stems from a known or unknown dynamics source. By leveraging a continual Bayesian Gaussian Mixture Model over a memory buffer, our framework can infer the identifiers and relationships of data over time - required for effective meta-learning. Empirical results on synthetic cardiac data demonstrate superior simulation forecasting, computational scalability, and resilience to catastrophic forgetting compared to existing baselines. Keywords: Cardiac EP, Personalization, Continual Meta-Learning

1. Introduction Personalized virtual models of cardiac electrophysiology, tailored to the observational data of individual subjects, have demonstrated utility in clinical applications such as treatment planning (Prakosa et al., 2018) and arrhythmia risk stratification (Arevalo et al., 2016; Trayanova et al., 2024). However, broader clinical integration faces two key challenges. First, calibrating these models to match an individual’s physiology remains a complex and time-intensive task, particularly for parameters that cannot be directly observed, such as tissue material properties (Niederer et al., 2020; Trayanova et al., 2024). Second, the high computational demands of running these simulations limit their scalability, making it difficult to embed them into clinical workflows or perform thorough uncertainty analysis (Niederer et al., 2020; Cluitmans et al., 2024). There has been considerable advancement in personalizing cardiac model parameters (Miller et al., 2021; Dhamala et al., 2018a; Wong et al., 2015; Sermesant et al., 2012; Zettinig et al., 2013; Cluitmans et al., 2024), particularly in the context of electrophysiological (EP) modeling. Early approaches typically involved iterative optimization schemes that required repeated evaluations of the computationally intensive EP simulation model (Wong et al., 2015; Sermesant et al., 2012). More recent efforts have incorporated modern machine learning (ML) techniques, including active learning (Dhamala et al., 2020), reinforcement learning (Neumann and Mansi, 2020), and direct ML of the parameter-to-observation mapping (Coveney et al., 2021; Giffard-Roisin et al., 2016). While these approaches have reduced the number of simulation calls required during optimization, the overall process of personalization remains complex and resource-intensive. Furthermore, even after obtaining a personalized model, the computational burden of running the simulation remains a significant obstacle for clinical-scale deployment. Simultaneously, the field of deep learning (DL) has seen growing interest in building efficient neural approximations for computationally expensive scientific simulations (Kasim et al., 2021). However, progress on neural surrogate modeling for cardiac EP remains limited. Early results have primarily focused on 2D domains (Cantwell et al., 2019; Kashtanova et al., 2021) and 3D simulations of the left atrium (Fresca et al., 2021), with more recent work extending to full-field simulations of the both ventricles (Martinez et al., 2025) or whole-heart electromechanical simulations (Salvador et al., 2024). A central challenge in this space is the dependence of simulations on physio2

logical parameters, denoted θ, such as material properties. Most approaches aim to learn a single neural function f (θ) to approximate the behavior of a simulation model M(θ). This introduces two main issues. First, learning an accurate f (θ) requires a large dataset of paired samples {θi , M(θi )} simulated across the input space of θ, which is computationally challenging to generate. While the emergence of physics-informed neural networks (PINNs) alleviates this need of data by supervising the neural networks with the analytical partial differential equations (PDEs) of M(θ), it requires a distinct PINN to be trained from scratch for each particular value θi (and boundary/initial conditions for that matter) (Herrero Martin et al., 2022; Gomez et al., 2025). In general, in either approach, there is the difficulty of ensuring sufficient coverage of the parameter space during training, such that the trained surrogate is robust to out-of-distribution parameter values at practical deployment. Second, even if f (θ) is well-approximated, its application depends on having accurate patient-specific values of θ, which are often unknown or unobservable, in order to become personalized (Trayanova et al., 2024). Recently, a learning-to-learn approach has been proposed to address the challenges above associated with learning how to personalize a neural surrogate itself (Jiang et al., 2022). The guiding principle is being interested in a set of personalized neural functions, rather than one single generic neural function, as the simulation surrogate: therefore, instead of learning a single neural surrogate, one can learn the process of learning a personalized neural surrogate from a small number of available data of a subject (context data). This is cast as a novel formulation of few-shot generative modeling via Bayesian meta-learning (Jiang et al., 2022). It includes two elements: 1) a set-conditioned generative model as the neural surrogate for cardiac EP simulation that, conditioned on an abstract latent embedding that represents patient-specific information, learns to generate target simulations that are personalized to an individual, and 2) a meta-model of amortized variational inference (VI) that learns to extract such patient-specific embedding via feed-forward embedding from context data of variable sizes from a patient. Compared to optimization-based meta-learning methods (Finn et al., 2017a; Ravi and Larochelle, 2016), this type of feed-forward meta-models remove the need of further training or fine-tuning, obtaining a model at meta-test time via simple feed-forward embedding of context data. With this, at test time, a personalized neural surrogate can be quickly obtained by simple feedforward embedding of a small and flexible number of data available from a 3

subject. Despite significant advancements in neural surrogates and their personalization - particularly through meta-learning approaches - these methods share a fundamental limitation: the assumption of access to a static training distribution of diverse subjects. When fast adaptation to a subject within or near the training distribution is possible, if data from novel subjects becomes available, such models require expensive retraining with the entirety of previous data to avoid catastrophic forgetting - the well-known phenomenon where models forget earlier tasks when trained only on new information (Kirkpatrick et al., 2017). In clinical practice, subject data typically arrives sequentially and often exhibits distributional shifts, rendering models trained on prior cohorts insufficient for generalization to new subjects or even known subjects changing over time (Kolk et al., 2024). Consequently, there is a critical need within the clinical workflow for efficient simulation surrogates that can be rapidly personalized to individual subjects while maintaining performance across non-stationary distributions of cardiac data (Bhagirath et al., 2024). The concept of continual meta-learning (CML), which combines continual training with few-shot adaptation, is well-suited for addressing this challenge. Existing CML methods, however, face two core limitations in this context. Most works rely on gradient-based meta-learners like model-agnostic metalearning (MAML) (Finn et al., 2017b), which leverages a shared parameter initialization to adapt to context data in a few gradient-update steps. This requires computationally expensive test-time fine-tuning and can fail to generalize beyond image classification and low-dimensional regression tasks (Mishra et al., 2017). Second, most employ continual strategies that are task-agnostic, leveraging the approximate equivalence between continual and meta-learning objectives in aligning gradients by directly applying the metamodel adapted on actively-streaming data to all prior tasks (Riemer et al., 2018; Joseph and Gu, 2021). These approaches fundamentally conflict with the requirements of mesh-based cardiac data, where per-vertex spatial variations induce patient-specific dynamics (that a shared initialization cannot accommodate), and where clinical deployment demands low-latency inference. Addressing these gaps requires moving beyond gradient-based and task-agnostic CML methods toward approaches that infer task identities directly and can adapt without expensive gradient updates. In this paper, we present a novel concept to achieve continually-adapted personalized neural surrogates in a single coherent end-to-end framework of 4

continual meta-learning (CoMetaPNS), addressing the challenge of continual adaptation central to personalized cardiac modeling. Specifically, we develop a continual Bayesian meta-learning framework that (1) learns to learn a personalized neural surrogate from limited, variable-sized subject-specific observations; (2) incorporates a task-relational past-sample reservoir that uses continual Bayesian Gaussian Mixture Models to distinguish whether new data stems from previously encountered or novel dynamics; and (3) continually updates the meta-inference model to preserve past knowledge while integrating new information, thereby avoiding catastrophic forgetting. By enabling automatic task identification and directly inferring per-vertex conditioning variables, CoMetaPNS aligns with the spatial variability inherent in patient-specific cardiac meshes, advancing beyond the limitations of existing gradient-based and task-agnostic CML approaches and enabling continual and rapid personalization even under sequential and non-stationary data streams typical of clinical deployment scenarios. We evaluated CoMetaPNS on synthetic non-stationary cardiac data streams and its generalization to real cardiac data, highlighting the joint necessity of its meta-learning and continual learning components. We included comparisons to the following baselines: (1) individually-personalized cardiac simulation models optimized using established techniques (Dhamala et al., 2018b), 2) generic neural surrogates lacking set conditioning or meta-inference, and (3) alternative meta-inference approaches representing standard practices in continual meta-learning (CML) literature. CoMetaPNS achieved superior personalization and predictive error at significantly reduced computational cost compared to both conventionally-optimized simulations and existing meta-inference baselines. Importantly, we found that neither metalearning nor continual learning alone sufficed under non-stationary data conditions - only their integration provided stable performance. Additionally, we demonstrated the utility of the feed-forward embedding alongside a continual Bayesian Gaussian mixture model, which reliably distinguished between re-emerged (i.e., known) and novel data sources. This enabled control over updating the meta-model’s parameters or applying fast adaptation alone. 2. Problem Formulation Consider a cardiac electrophysiology simulation model x1:T = M(v; θ), where v is a known input (e.g., electrical stimulation applied to the virtual heart), θ is an unknown parameter representing subject-specific properties 5

(e.g., tissue conductivity), and x1:T are heart surface unipolar potential maps. The objective of model personalization is to estimate θ̂ that minimizes the discrepancy between observed and simulated x1:T . A personalized neural surrogate aims to learn a function f (v, θ̂) that both 1) approximates the simulation output M(v; θ̂) for computational efficiency and 2) aligns with the individual’s observations for personalization. A natural approach is to learn a neural function f (v, θ) that approximates the simulation model M(v; θ) across the space of θ (Kashtanova et al., 2021; Fresca et al., 2021), yielding a generic surrogate that either requires an explicit input of θ or a separate optimization step for personalization. In contrast, leveraging state-of-the-art meta-learning, we aim to derive a set of personalized neural surrogates that automatically adapt to a small, variableκ sized set of context observations X s = {xs,i 1:T }i=1 from an individual subject, where κ is the number of observations available for the subject. This is formulated in the framework of Bayesian meta-learning as: Z s p(x̂1:T |v, X ) = p(x̂1:T |v, c) qζ (c|X s ) dc, where p(x̂1:T |v, X s ) describes the likelihood of the sequence x̂1:T given known input v and context observation sequences X s . The generative model p(x̂1:T |v, c) conditions on known stimulation input v and a latent embedding c personalized to the subject, while qζ (c|X s ) is a meta-model parameterized by ζ that learns to extract patient-specific c through feed-forward embedding of the context observations X s . We aim to maintain the performance of the meta-model over a nonstationary stream of observations as data from different subjects become available over time, which consists of varying heart meshes and material properties θ. To reflect practical constraints in clinical settings, we make the following assumptions regarding the nature and availability of such data: 1. We do not assume access to the complete historical data distribution as it is often infeasible due to memory/storage limitations and privacy concerns. Instead, we assume access to a fixed-size memory buffer that can hold up to M samples. 2. We assume that we know when the data source - i.e., the subject changes, a reasonable assumption within clinical workflows. We however do not assume access to ground-truth identifiers for these sources: thus, we do not know the similarity or dissimilarity about the underlying conditions among subjects; similarly, we do not know whether 6

Figure 1: Overview of CoMetaPNS, showing A) the framework of few-shot generative modeling via Bayesian meta-learning that B) continually aggregate a heterogeneous data stream of cardiac electrophysiology dynamics with a sample reservoir that identifies subjects via Gaussian mixture models.

any newly-presented data corresponds to a previously-known or novel subject. 3. We assume the presence of local stationarity where data from each subject are presented for some period of time. This is a common assumption utilized in continual-learning literature to allow the neural model to optimize sufficiently and stably before a new data source has the potential of appearing. Within these constraints, our goal is to effectively manage the limited memory buffer to approximate the overall data distribution over time and infer the relationship between incoming and previously encountered data sources, such that we can learn to adapt our neural surrogates to data from such non-stationary distributions of subjects without forgetting. 3. Methodology Fig. 1 gives an overview of CoMetaPNS. It includes two major components: A) an adaptive neural surrogate realized via feed-forward metalearning, and B) a continual learning strategy for learning this meta-model across continuous data streams of diverse cardiac EP dynamics. 3.1. Adaptive Neural Surrogates via Feedforward Meta-Learning The adaptive neural surrogate is similar to that described in (Jiang et al., 2022), which includes a set-conditioned generative model (Section 3.1.1) and 7

a meta-inference strategy for extract patient-specific embedding to condition this generative model (Section 3.1.2). 3.1.1. Set-Conditioned Generative Model The generative neural surrogate is conditioned on a context-set embedding c and known stimulation input v. It consists of two components: a temporal transition model Tθt that governs the evolution of the latent state zt , and a spatial emission model Gθs that maps the latent state to the highdimensional cardiac mesh output xt . These components are defined as: Transition: zt+1 = Tθt (zt ), Emission: x̂t+1 = Gθs (zt+1 ). The initial latent state z0 is derived from the known stimulation v using a neural function composed of a linear layer applied to the embedding of v, denoted as z0 = fρ (v). Temporal transition Tθt : The temporal transition function Tθt is inspired by Gated Recurrent Units (GRUs) (Chung et al., 2014). Note that, with regular GRUs, Tθt would be global to the training data rather than subject-specific. Instead, we condition Tθt on the context-set embedding by creating conditional gated transition functions: (1)

= ELU(α1 zt + β1 c + γ1 ),

(2)

= ELU(α2 zt + β2 c + γ2 ),

zt

zt

(3) zt = α3 zt + β3 c + γ3 , (1) gt = σ(W1 zt + b1 ),

(1)

(2)

ht = ELU(W2 zt + b2 ), (3)

zt+1 = (1 − gt ) ⊙ (W3 zt + b3 ) + gt ⊙ ht , where c ∼ p(c|X s ) (see Section 3.1.2) and θt = {Wi , bi , αi , βi , γi }3i=1 are learnable parameters. The model has flexibility to choose a linear transition for some dimensions and non-linear transition for the others. Spatial emission Gθs : Since x resides on the 3D geometry of the heart, we implement the emission function Gθs using Graph Convolutional Neural Networks (GCNNs). We represent the heart’s triangular mesh as an undirected graph, where edges between vertices are assigned attributes defined by the normalized differences in their 3D coordinates, provided an edge exists. 8

Encoding and decoding are performed over hierarchical representations of this geometry, constructed via a specialized mesh coarsening algorithm (Cacciola, 2007). To enable spatial convolution across graphs, we adopt a continuous spline kernel (Fey et al., 2018). To enhance the model’s expressivity and depth, we incorporate residual blocks using skip connections implemented with 1D convolutions, following (Jiang et al., 2020). 3.1.2. Meta-Model for Amortized Variational Inference We introduce a meta-model pϕ (c|X s ) to model the conditional distribution of c given a context set X s , realized via a feed-forward neural network hϕ . First, each sample xs1:T ∈ X s is embedded through the neural function hϕ (xs1:T ) that uses a GCN-GRU cell (Jiang et al., 2021) to obtain the sequential information from the graph, and aggregate it across time with a linear layer. We then average all latent embedding in Xks : 1 X hϕ (xs1:T ), s ∈X s s x |X | 1:T

(2)

which then parameterizes pϕ = N (µc , σc2 ) via two separate linear layers. The conditional factor c is then sampled by c = µc + ϵ ⊙ σc , where ϵ ∼ N (0, I) (Kingma and Welling, 2013). Now consider a set Q of subjects {Qk }K k=1 , where the data distribution of each subject Qk is defined by p(Qk ) and the distribution of subjects by p(Q). For each subject Qk , we define two associated sets of high-dimensional q,2 q,Nk cardiac data: 1) the query set Xkq = {xq,1 1:T , x1:T , ..., x1:T }, for which only the stimulation inputs Vk = {v 1 , v 2 , ..., v Nk } are available to forecast from; and s,2 s,Mk 2) the context data Xks = {xs,1 1:T , x1:T , ..., x1:T }, consisting of fully-observed sequences used to infer the subject-specific conditioning variable c. Note that Mk ≪ Nk and, as discussed below, these sets are mutable during episodic training. The evidence lower bound (ELBO) we optimize for each Qk can be for-

9

mulated as: Nk X q,i k,i s I(xq,i 1:T , Qk ) log p(x̂1:T |v , Xk ) ≥ i=1 Nk X

h q,i k k,i I(xq,i , Q ) k Eqζ (ck |X s ∪xq,i ) [log p(x̂1:T |c , v )] 1:T k

1:T

i=1

i k s − KL(qϕ (ck |Xks ∪ xq,i )||p (c |X )) ϕ k 1:T (3) where I(x1:T , Qk ) a binary membership function that equal to 1 if x1:T is in Qk and 0 otherwise. Note that we let pϕ (ck |Xks ) and qϕ (ck |Xks ∪ xq1:T ) share the same meta set-embedding networks to parameterize their means and variances. We further regularize pθ (ck |Xks ) to be close to a standard Gaussian distribution N (0, I), giving a loss function LQk (Ik ) per Qk as: LQk (Ik ) =

Nk X

h   k k,i Ik Eqϕ (ck |X s ∪xq,i ) log p(x̂q,i 1:T |c , v ) k

1:T

i=1

k s − λ1 KL qϕ (ck |Xks ∪ xq,i 1:T )||pϕ (c |Xk )
i − λ2 KL pϕ (ck |Xks )||N (0, I)

(4)

where Ik is shorthand for I(xq,i 1:T , Qk ), and λ1 and λ2 are regularization multipliers. Let Θ = {θs , θt , ρ, ϕ} including all unknown weight parameters of our model described above, we have the overall optimization objective across T as: X Θ̂ = max LQk (Ik ) (5) Q ∈Q,Q∼p(Q) Θ

k

The optimization problem in Equation 6 is solved in episodic training. In each training episode across all subjects, the input data is divided into two separate sets: a context set Xks consists of small sets of samples from each subject and the query set X q formed by the remaining data. The model is asked to take Xks for each subject k, derive ck , and generate samples X̂kq from the target set stimulation inputs Vk . 3.1.3. Use of Feed-Forward Adaptation A natural question arises regarding the use of a feed-forward meta-learner to perform the adaptation within this setting rather than standard gradientbased meta-learners. The most common realization in CML is via the modelagnostic meta-learner (MAML) (Finn et al., 2017b), in which the meta-model 10

is adapted to new tasks by optimizing a learned parameter initialization using a few gradient steps on the available context data. While effective in classification tasks, gradient-based meta-learners like MAML often struggle to adapt to the optimization landscapes found in other domains (Mishra et al., 2017). Moreover, in the continual environment, they require fine-tuning for every new task, including previously encountered ones, making them sensitive to optimization hyperparameters and computationally inefficient. This inefficiency becomes particularly problematic in cardiac simulation scenarios, where the underlying neural surrogate models are often inherently computationally expensive to update via backpropagation. Furthermore, the algorithmic prior of gradient-based meta-learners - the mechanism through which they adapt to task-specific information - conflicts with the demands of data represented as heterogeneous graph structures. In our feed-forward conditioning framework, per-vertex conditioning vectors can be generated in a way that remains agnostic to the global graph topology while still providing spatially localized adaptation when passed through a shared dynamics function. In contrast, gradient-based methods require a shared initialization across all tasks, adapted via a small number of gradient steps. In graph-based domains, where the structure varies between tasks, this shared initialization cannot capture per-vertex nuances, since vertex-specific parameters are not transferable across different graph topologies. Consequently, adaptation is restricted to the shared parameters of the dynamics function, making task-specific adaptation substantially more challenging. Within just a few gradient updates, the meta-learner must simultaneously adapt the shared dynamics to both the temporal behavior and the diverse spatial configurations of each new graph - an inherently difficult objective given the complexity of the domain. 3.2. Continual Learning: How to Adapt Without Forgetting Now consider a condition where the set of subjects Qk and their corresponding data are presented sequentially over time. Further consider that a subject’s data may reappear multiple times in this stream of data. To be general, we assume that we know when then source of the data (i.e., subjects) changes, but we do not know the identity of the source (i.e., we do not know if a newly-presented subject corresponds to a new or previously seen subject). In continual learning, this corresponds to a setting where taskboundary (i.e., the boundary of data-distribution switching) is known, but the identifier of the underlying task is unknown. 11

The learning objective as defined in Equations 6 is optimized over a stationary distribution of p(Q), as in standard meta-learning practice. As such, despite its ability to infer context and personalize the resulting neural surrogate, its meta-components become susceptible to catastrophic forgetting in non-stationary distributions of p(Q). Additionally, how to properly identify the underlying task to enable context-query data pairing within the task becomes complicated. Below, we describe mechanisms to approximate the stationary subject distribution p(Q) and the subject membership function Ik in order to achieve a continual approximation of the learning objective in Equation 6. These include a reservoir sampling based experience replay strategy (Section 3.2.1), combined with two alternative mechanisms of continual meta-learning depending on how task identifiers are estimated and their relations modeled. (Sections 3.2.2 – 3.2.3). 3.2.1. Experience Replay based on Reservoir Sampling In the non-stationary setting where only the data of one distribution is actively streaming in at a given time, a first question is how to obtain and aggregate errors from prior data distributions to update the meta-weights. For this, we adopt the reservoir sampling method commonly used in existing CML approaches, where a simple algorithm tracks the number of samples (N ) seen and, for each incoming sample, overwrites an existing buffer sample with probability M/N where M is the size of the reservoir. With actively streaming data and replayed samples in the reservoir buffer, the next question is how to accurately pair context and query samples from the same task in order to approximate the meta-objectives as described in Equation (6). Under the assumption of local stationarity, this question is trivial for the streaming data which is assumed to be the current task: as adopted in prior continual meta-learning (CML) works (Caccia et al., 2020; He et al., 2019), the most recent k observations are used as the context set for the current task. For samples in the reservoir buffer, we consider two mechanisms for estimating their task identifiers and model the relations among the identified tasks. Depending on the mechanisms used to estimate the task identifiers, we also adjust the reservoir samples to approximate p(Q) differently which will be detailed in the next two sections. 3.2.2. Task-Aware Meta-Learning A simple mechanism is to directly leverage the known boundaries of subject switching, and assume that the incoming data source at every boundary 12

belongs to a unique novel task. Under this assumption, we can simply use a boundary counter to assign pseudo task labels to data samples. When sampling from the reservoir, context-query pairs are easily matched by using each sample’s pseudo task label. Because we know the boundary of subject switching but not the identifiers of the subjects, in relation to the original meta-learning objective in Equation (6), this approach may represent a subject Qk with one or multiple tasks, where each task will correspond to at k most one subject, i.e., Qk = {Tk,i }m i=1 , where mk represents the number of times the subject reappeared in the data stream. The membership function Ik,i = I(x1:T , Tk,i ) for a sample x1:T in task Tk,i is exact given the known task boundary. In other words, in this approach, the meta-loss LQk (Ik ) per subject Qk in Equation (6) will be approximated as: Xmk LTk,i (Ik,i ) LQk (Ik ) ≃ (6) i=1

and the overall meta-loss over the streaming and replayed samples becomes: X LQk (Ik ) (7) Q ∈R,R∼p(R) k

where R includes limited samples of each Qk in the reservoir, and p(R) approximates p(Q) in the reservoir. This approach has the advantage of being straightforward and easy to implement. Furthermore, there is no error in estimating the task membership of a sample, which will avoid pairing context and query samples from different data distributions. However, by overlooking the underlying relationships among tasks, this simplified approach is likely to consider data from similar and even re-appeared subjects as distinct tasks. In clinical settings, this can result in unnecessary retraining on familiar subjects instead of utilizing existing meta-knowledge, leading to wasted time and computational effort. Furthermore, as the number of subject-switching boundary increases continually, this method can result in an increasingly-fragmented reservoir buffer. Given the finite size of the reservoir, each task will be represented by a reducing number of samples: this can result in overfit to limited contextquery pairs per task and, when the data distribution of some subjects are very different from the rest, they will become under-represented tasks in the reservoir and risk being rapidly forgotten. Collectively, this will result in an increasingly inaccurate p(R) as an approximation of p(Q) in Equation (8) as the number of tasks increase with the number subject-switching boundaries. 13

3.2.3. Task-Relational Meta-Learning An alternative mechanism is to explicitly model task relations based on similarities among the streaming and replayed data samples. To this end, we leverage a continual Bayesian Gaussian mixture model (GMM) to describe the context-embeddings from the reservoir samples as a distribution of clusters p(R), where each cluster is defined as Rk , within the reservoir over time and determines a sample’s relation to their data distributions p(Rk ).The GMM is updated continually by re-initializing its parameters using the cluster components and weights from the previous task, ensuring temporal consistency in clustering. This approach removes the dependence on traditional reservoir sampling strategies for memory updates; instead, the GMM itself determines which samples to retain or overwrite based on inferred task similarity. The continual GMM integration proceeds in four stages, assuming we are at a newly encountered boundary at the j-th subject switch. At this point, an initial batch of new, unlabeled samples [Xkq , Xks ] is observed, and the GMM-based procedure for clustering and memory update begins as follows: i) Using the previously fit GMM, we compute the embeddings csk of the current batch’s context set Xks and evaluate their per-cluster log-likelihoods. This provides a quantitative measure of whether the new data distribution aligns with previously observed sources. If the average log-likelihood of these samples is negative, we classify the data source as novel and proceed by updating the meta-model’s parameters through full backpropagation. Novel data sources are assigned a unique task ID, incremented based on the number of unique tasks identified so far. Conversely, if the average log-likelihood is positive, we classify the source as known and opt for fast feed-forward adaptation without parameter updates. To assign the task ID for this known data source, we identify the cluster to which its samples belong, then assign it the majority task ID from the reservoir Cjreservoir ’s samples in that cluster. Note that we do not use the clusters themselves as task identifiers, as this could lead to per-sample misclassifications if a sample from a given data source is incorrectly clustered. Instead, we aim to maintain stable task ID assignments over time by grouping data sources consistently, while still leveraging the GMM’s cluster likelihoods to perform effective novelty classification. Overall, this identification mechanism enables dynamic control over when to apply efficient adaptation versus full meta-model training, based on inferred novelty, as further demonstrated in Section 4.5.1. 14

ii) During the course of task j, a subset of its streaming samples is stored in an auxiliary memory buffer Cjactive using the standard reservoir sampling procedure. This buffer serves as a temporary holding set and is maintained independently of whether the current task has been identified as novel or known. The size of Cjactive is strictly smaller than that of the main reservoir and can be dynamically adjusted as a fraction of the current number of 1 , to ensure balanced representation as the number of clusters, specifically ||R|| identified sources grows. iii) When the (j + 1)-th subject-switching boundary is encountered, prior to processing any new incoming data, the GMM is re-fitted using the context embeddings from both the existing reservoir samples Cjreservoir and the auxiliary buffer Cjactive collected during step ii. This re-fitting occurs regardless of the task classification outcome in step i (novel or known). Specifically, a new mixture component is initialized from the mean and covariance of Cjactive , and assigned a low initial weight. Upon re-fitting, the updated weight of this component reflects the novelty of the data: high weight indicates a novel source not yet represented in the reservoir, while low weight suggests redundancy with existing clusters. To prevent unbounded growth of the mixture model, components with weights below 0.05 are pruned. iv) After fitting the GMM, but still prior to new data processing, the next step is to integrate the auxiliary buffer’s samples into Cjreservoir while ensuring two objectives: 1) preserving the reservoir buffer size and 2) maintaining overall task distributional balance. This integration strategy is guided by the novelty classification from step i. If the task is classified as novel, we sample an equal number of entries per unique task from the union of the reservoir Cjreservoir and auxiliary reservoir Cjactive , ensuring balanced representation across all tasks. If the task is identified as known, we restrict the update to rebalancing only the corresponding task’s entries within the reservoir, allowing new samples to be integrated without disrupting the broader distributional structure. Finally, in relation to the original meta-learning objective in Equation (4), this approaches treats each GMM cluster Rk as a unique task, assigning each cluster its own pseudo-label for its assigned reservoir samples that spans across subject-switching boundaries. Given this is a data-driven clustering process, it is possible that a cluster may consist of data samples stemming from more than one underlying subject. This is a beneficial aspect of the clustering in which subjects whose underlying data distributions overlaps significantly are mapped to a similar context-embedding space. This allows 15

the meta-model to consider the relationship among subjects and prevent unnecessary retraining on similar data distributions. Formally, the data distribution as described by the GMM within the reservoir buffer, including the distribution p(Rk ) of each cluster Rk and the distribution of the clusters p(R), approximate the actual data distribution of each subject p(Qk ) and subject distribution p(Q), respectively. This gives a continual approximation of the meta-objective in Equation (6) as: X LRk (Mk ) (8) R ∈R,R∼p(R) k

where LRk (Mk ) is calculated over the distribution of p(Rk ) and the estimated cluster assignment Mk approximates the true subject membership Ik . 4. Synthetic-Data Experiments In all experiments, CoMetaPNS consisted of 4 GCNN blocks and 2 regular convolution layers in the encoders, 1 context-set aggregator with a GCN-GRU block followed by a linear layer to compress time and another linear layer for feature extraction, 1 conditional gated transition unit for the generation of latent dynamics, and 4 GCNN blocks and 2 regular convolution layers in the decoder. We used Adam optimizer (Kingma and Ba, 2014). The learning rate at the beginning of each task was set at 1 × 10−3 with a learning rate decreasing at a rate of 0.5 every 100 iterations. The two KL multipliers were: λ1 = 10−4 and λ2 = 0.1. All experiments were run on Tesla T4s with 16GB memory. Our implementation is available upon acceptance. 4.1. Experimental Setup 4.1.1. Data In synthetic data, we considered 3 heart meshes comprising 448, 475, and 480 nodes, respectively, with a combination of 12 distinct spatial configurations of scar tissue representing various injury patterns. The anatomical information for these heart meshes were derived from real subject hearts, sourced from Dawoud (2009) and anonymized here. Each configuration was treated as a separate subject. For the task sequence, we selected 9 subjects (taking 3 scars from each heart) in random order, yielding a total of 517 unique simulations. The remaining three scar configurations, totaling 41 simulations, were held out to evaluate generalization performance following continual learning. 16

On each subject, we simulated macroscopic spatiotemporal propagation sequences of action-potentials by the Aliev-Panfilov model (Aliev and Panfilov, 1996) considering approximately 57 ± 11 different origins of activation separated into a 80-20 train-test split. To obtain the heart-surface measurements in the form of extracellular potential from volumetric action potential, we obtained the forward operator by solving Poisson’s equation using the coupled mesh-free method and boundary element methods as described in (Wang et al., 2009; Liu, 2003; Brebbia et al., 2012). Because the action potential simulated by the Aliev-Panfilov model was unit-less in both amplitude and time, the generated signals on the heart were also unit-less. Specifically, in our experiments, we considered the depolarization process that was downsampled in time and represented by 55 discrete time steps. The synthetic simulation data for the heart meshes is made available upon acceptance. 4.1.2. Baselines To the best of our knowledge, there is no existing work that considers continual learning of cardiac simulation models. As such, there are no direct baselines to compare against. Instead, we adapted relevant stationary models as best as possible to evaluate the contribution of CoMetaPNS. Non-learning-based personalization methods. We included classic nonlearning baselines that apply iterative optimization routines to estimate the parameters of a cardiac simulation model. In specific, we considered the published work where tissue excitability within the Aliev-Panfilov model was estimated by derivative-free Bayesian optimization (Dhamala et al., 2018b). We included both optimization formulations discussed in (Dhamala et al., 2018b): one parameterizing tissue excitability using a seven-segment division of the cardiac mesh (FS-BO), and another leveraging a variational autoencoder (VAE) to parameterize spatially-varying tissue excitability (VAE-BO). Since these methods were designed to optimize a single subject-specific simulation model at a time, for each of the 9 subjects, we optimized the AlievPanfilov model using the same context data used to adapt CoMetaPNS. The resulting Aliev-Panfilov model was used to simulate/predict the action potential sequence using other origins of activations in the query set. Learning-based neural Surrogates. While neural surrogates for cardiac simulations have been explored (Cantwell et al., 2019; Kashtanova et al., 2021; Fresca et al., 2021), existing approaches are typically limited to 2D models on image grids or focus on atrial structures (Fresca et al., 2021). 17

These models, once trained, remain general and require additional optimization to personalize predictions for subject-specific data. To more directly evaluate the contributions of our meta-model against traditional neural surrogates, we implemented a version of CoMetaPNS without the meta-model or continual learning - resulting in a standard generative model p(x1:T |s, c), where the embedding c is inferred directly from y1:T via q(c|y1:T ). We refer to this baseline as PNS. This comparison allows us to isolate and assess the impact of the meta-learning and continual learning components individually. Alternative CML baselines. To benchmark against traditional CML approaches, we include a gradient-based variant of our meta-model, denoted as MAML-PNS. This model adopts the standard MAML framework, where the shared initialization parameters correspond to the learnable weights of the GRU-based transition function T . Specifically, it optimizes over the parameter set {Wi∗ , b∗i , αi∗ , βi∗ , γi∗ }3i=1 , which are updated through a small number of gradient descent steps using the available context data for each task. 4.1.3. Evaluation Settings & Metrics To assess the effect of non-stationary subject distributions on CoMetaPNS and its baselines, we evaluated PNS, MAML-PNS, and CoMetaPNS under three experience replay settings: i) Naive Learning (NL), where no past samples are replayed; ii) Exact Replay (ER), in which the reservoir can accommodate all past samples; and iii ) Task-Aware Meta-Learning, as introduced in Section 3.2.2. For CoMetaPNS, we further delved into a comparison between Task-Relation Meta-Learning and Task-Aware Meta-Learning strategies. Since the non-learning FS-BO and VAE-BO baselines performed optimization for each subject instance independently, they were not affected by this setting of data streams. The accuracy of the solutions in forecasting was measured by the mean square error (MSE), spatial correlation coefficient (SCC), and temporal correlation coefficient (TCC) between the reconstructed and actual potential sequence on the heart surface. In addition, we considered the dice coefficient (DC) of the abnormal tissue region obtained by thresholding signals with Otsu’s method (Otsu, 1979). To evaluate the ability of a model to learn continually, we followed (Riemer et al., 2018) and examined the above forecasting metrics in two perspectives: Retained Performance (RP) and Learning Performance (LP). The RP metric represented the average performance across all tasks after they were sequen18

tially considered, emphasizing a model’s ability to retain performance on previous tasks. The LP metrics reflected the average performance on a task immediately after it was learned, measuring a model’s effectiveness in incorporating new information. To assess the extent of catastrophic forgetting, we also reported the average difference between LP and RP metrics, referred to as Backward Transfer and Interference (BTI), with negative values indicating increased forgetting (Riemer et al., 2018). When comparing the choice of the meta-learner in Section 4.3.3, we included two additional metrics related to computational efficiency. Time-toAdapt-N (TTA-N) represented the average time it takes, in seconds, for a model to adapt to its context data in the presence of N tasks. Time-to-Train (TTT) represented the average time it takes, in minutes, for a model to train over a sequence of tasks. 4.1.4. Other Implementation Details All experiments were run on NVIDIA Tesla T4 GPUs with 16GB memory in instanced cloud systems to control hardware purity. We used PyTorch 1.13.1 and scikit-learn 1.4.2 for deep learning optimization and GMM fitting, respectively. All models were trained to forecast the 55 timesteps using only the onehot mesh encoding of the activation origin. We set the window of local stationarity to 2000 iterations, where each iteration is composed of 5 active task samples. For CL, we used a reservoir size that can accommodate approximately 45% of the total data at any given time, resulting in a reservoir size of 200. The size of the k-shot context set was set to vary, where at each iteration we sample k between 1-5. To ensure fairness, comparable model components between the meta-learners were scaled to maintain consistent total parameter counts of approximately 1.1 million trainable parameters. Shared backbone components had identical hyper-parameters, while modelspecific hyper-parameters were tuned over the same amount of GPU hours. To test LP and RP metrics, we considered held-out testing sets of every dynamics. To test LP, whenever a subject was finished training for the first time, it was evaluated on its testing set for its performance. At the end of the task sequence, on the resulting model, we evaluated it on every subject independently and averaged their performances to get the RP metric. BTI, then, was the average difference between all subjects’ LP and their respective RP.

19

4.2. The Effect of Continual Learning In this section, we examine the necessity of employing suitable continual learning strategies to address the challenges posed by non-stationary cardiac dynamics. 4.2.1. Quantitative Results We present results on spatial correlation coefficient (SCC) for the synthetic heart dataset in Fig. 2. Gray and blue shaded regions correspond to worst- and best-case scenarios, respectively, representing models trained without (NL) or with full (ER) access to prior task examples. Red (VAEBO) and green (FS-BO) dotted lines mark the average performance of classic optimization methods personalized to each heart, serving as a reference for RP performance. The blue dotted line represents the meta-model trained in a full stationary setting with access to the full data distribution, termed MetaPNS. For LP, all models showed relatively stable performance across experience replay settings, suggesting minimal impact of replay on learning the active task. Interestingly, CoMetaPNS exhibited slightly better LP under the naive setting, potentially due to fully utilizing model capacity on the current subject without the influence of replay. In contrast, RP results revealed that all models suffered from performance degradation when replay was absent, indicating severe forgetting. Introducing a memory mechanism - whether through exact replay or task-aware sampling - consistently mitigated this issue and, for CoMetaPNS, even display positive forward transfer on early tasks. Among all, the proposed combination of a meta-model and task-aware meta-learning provided the most substantial gains, highlighting the necessity of integrating both continual and meta-learning components. Fig. 3 further illustrates per-task performance for CoMetaPNS under settings with (left) and without (right) task-aware meta-learning strategies. Overall, we can see that CoMetaPNS consistently outperformed all baselines considered including the classic optimization approaches, and matched the proposed meta-model trained in a stationary setting. 4.2.2. Visual Results Fig. 4 provides visual examples of the forecasted electrical activity from the proposed meta-model under each continual strategy, compared against the ground truth. After completing training on the final task, the metamodels were evaluated on the first task when it reappeared at meta-test time. 20

Figure 2: Spatial correlation coefficient (SCC) performance comparison on the synthetic data for the considered baselines across the continual strategies. Horizontal dotted lines represent the average performance of the per-subject baselines. All methods were run over 5 shared seeds.

Figure 3: CoMetaPNS with Task-Aware Meta-Learning (left) overcomes catastrophic forgetting in the continual setting vs. naive training (right).

Both reservoir sampling and exact-replay strategies produced accurate forecasting of the propagation dynamics, demonstrating effective mitigation of catastrophic forgetting. Notably, due to variability in subject data availability, exact-replay required sampling uniformly across task IDs rather than the reservoir itself to prevent degraded performance on underrepresented tasks. In contrast, the naive strategy failed to recover the subject-specific tissue properties, resulting in poor reconstruction performance. For the PNS baselines, we observe that while they effectively capture the initial propagation of activation, they fail to generate smooth or accurate dynamics in the latter half of depolarization. We attribute this to their ability to learn a generic representation of propagation dynamics, lacking the flexi21

Figure 4: Visual examples of reconstructed electrical activity from the proposed continual meta-model when using continual strategies. Both epicardial and endocardial surfaces are presented. The color bar shows the scaled range of the signal since it is unitless in synthetic data. The MSE value is shown for each model at each timestep. C: CoMetaPNS, M: MAML-PNS, P: PNS. ER: Exact-Replay, RS: Reservoir Sampling, NL: Naive Learning.

bility to adapt to subject-specific features. In contrast, MAML demonstrates more coherent long-term forecasting in terms of per-node activation patterns, but its predictions exhibit notably weaker activation strength. We attribute this to the inherent challenge of simultaneously adapting both spatial and temporal dynamics from a shared initialization in this domain. 4.3. The Effect of Meta-Learning In this section, we investigate 1) the necessity of the meta-model in learning across heterogeneous cardiac dynamics (by comparing with the PNS baseline without meta-learning), and 2) the impact of different algorithmic priors 22

Figure 5: Dice coefficient (DC) performance comparison on the synthetic data for the considered baselines across the continual strategies. Horizontal dotted lines represent the average performance of the per-subject baselines. All methods were run over 5 shared seeds.

used for meta-adaptation (by comparing with the MAML-PNS baseline that uses MAML as the meta-learner). 4.3.1. Quantitative Results Fig. 5 presents Dice coefficient (DC) results on the synthetic heart dataset, with full metric comparisons in Table 1. A clear performance gap emerges between non-meta and meta-learning approaches: the non-meta PNS baseline fails to consistently capture subject-specific features, while both metalearners demonstrate improved localization and forecasting. Comparing the two meta-learning strategies, feed-forward CoMetaPNS outperforms MAMLPNS, which only reaches parity with classical optimization methods. This distinction is particularly visible when comparing performance across metrics. Spatial correlation coefficient (SCC) captures overall signal trends over the mesh but is less sensitive to localized boundaries, such as those defining scar tissue. As a result, SCC may overstate the performance of models that align with global signal structure but lack spatial specificity. The Dice coefficient, in contrast, directly evaluates spatial overlap in activation maps and better reflects a model’s ability to recover subject-specific anatomical detail. Given these observations, we interpret the performance trends of MAMLPNS as indicating an ability to adapt to subject-specific characteristics and represent inactive tissue, but a limitation in capturing global propagation patterns. This aligns with our rationale in Section 4.1.2 regarding the difficulty of task adaptation imposed by the gradient-based algorithmic prior in this domain. Conversely, the generic PNS model captures global dynamics 23

Table 1: Performance metrics of 1) CoMetaPNS, 2) MetaPNS, 3) MAML-PNS, 4) PNS, 5) FS-BO, and 6) VAE-BO during meta-test. For the continual methods, metrics are evaluated on the model after continual training with the Task-Aware Meta-Learning strategy. Classical optimization methods remain fit per-subject as a baseline.

Model CoMetaPNS MetaPNS MAML-PNS PNS FS-BO VAE-BO

MSE

SCC

DC

4.3±1.0e-4 0.79±0.05 0.83±0.07 4.5±1.2e-4 0.74±0.09 0.86±0.05 6.7±0.9e-4 0.40±0.06 0.48±0.06 4.5±1.4e-4 0.49±0.06 0.28±0.05 5.3±6.5e-4 0.69±0.25 0.48±0.34 4.8±2.5e-4 0.46±0.15 0.48±0.09

effectively through shared parameters across subjects but lacks the capacity to localize subject-specific features due to its static nature. In contrast, the feed-forward conditioning approach demonstrates its advantage through pervertex conditioning, achieving both strong global dynamics alignment and accurate localization of subject-specific features such as scar regions. 4.3.2. Visual Results Fig. 6 presents visual comparisons of the forecasted electrical activity from the evaluated meta-learners - feed-forward and gradient-based - against the non-meta-learning baseline. Following training on the final task, each meta-model was assessed on the fifth task during meta-testing. The results demonstrate that meta-learning is essential for producing accurate personalized neural surrogates: the non-meta baseline (PNS) fails to recover subject-specific tissue properties and exhibits limited forecasting capability. Furthermore, the selection of the meta-learning algorithm significantly influences the capacity for personalization. The gradient-based MAML model shows weak activation forecasting throughout the sequence and does not sufficiently adapt to subject-specific characteristics in this instance. 4.3.3. Computational comparisons Table 2 compares the computational performance of the two meta-learners for model adaptation/personalization, versus classic instance-wise Bayesian optimization as a reference. As shown, the feed-forward adaptation adopted in MetaPNS substantially outperformed MAML in both the time it takes to 24

Figure 6: Visual examples of reconstructed electrical activity on both meta-learners considered alongside a non-meta baseline. Both epicardial and endocardial surfaces are presented. The color bar shows the scaled range of the signal since it is unitless in synthetic data. The MSE value is shown for each model at each timestep. C: CoMetaPNS, M: MAML-PNS, P: PNS. ER: Exact-Replay.

adapt a model at inference time (TTA) and the time it takes to train (TTT), across all continual strategies. Note that MAML exhibited a linear increase in adaptation time as the number of unique tasks grew. In comparison, the feed-forward approach maintained strong efficiency due to its parallelizable inference mechanism. When compared to classic optimization baselines, CoMetaPNS achieved similar or better predictive performance at a dramatically lower computational cost: a single Aliev-Panfilov simulation averaged 5 minutes, while the neural surrogate required only approximately 0.71 ± 0.05 seconds. Moreover, Bayesian optimization in the FS-BO required approximately 100 simulation calls per personalization, whereas CoMetaPNS completed the context set embedding in 0.49 ± 0.05 seconds on average. This computational advantage of the feed-forward meta-learner, especially its adaptation efficiency at inference time compared to both classic perinstance optimization and alternative MAML-type of meta-learners, make it appealing for clinical deployment. 4.4. Meta-Model Pre-Training In many clinical deployment scenarios, it is not uncommon to expect an initial stage of data collection and model training prior to model deployment. 25

Table 2: Complete adaptation efficiency comparison across baselines. Meta-learning methods were adapted over 100 batches of data across the synthetic heart dataset, timing each batch independently. Bayes-Opt represents the average per-subject time to personalize for the FS-BO baseline. TTA [X]: Time-to-adapt X tasks. TTT: Time-to-train.

Model

Metric

Time

TTA [1] (seconds) Bayes-Opt TTA [9] (seconds) TTT (minutes)

∼ 30000 ∼ 360000 ∼ 6000

TTA [1] (seconds) TTA [9] (seconds) TTT (minutes)

1.37 ± 0.08 4.17 ± 0.12 3789 ± 108

TTA [1] (seconds) Feed-Forward TTA [9] (seconds) TTT (minutes)

0.49 ± 0.05 1.48 ± 0.07 696 ± 23

MAML

This opens the possibility that continual learning does not begin entirely from scratch but rather from an initial distribution of subjects. Indeed, some CML works consider this scenario by first pre-training a meta-model on a small stationary subject distribution prior to its use in a continual setting (Caccia et al., 2020; Harrison et al., 2020). In this section, we considered CoMetaPNS in this setting and evaluated the effect of this meta-model pre-training on the sample efficiency of CoMetaPNS in continually adapting to new cardiac dynamics. More specifically, we considered three additional held-out scar-geometry configurations in addition to the nine considered in the previous sections. We started with a MetaPNS pretrained a stationary distribution of the nine tasks, and compared 1) its performance if directly applied to the held-out tasks (meta-generalization), 2) its continual learning over the held-out tasks (CoMetaPNS-pretrained), and 3) CoMetaPNS trained from scratch on the same tasks (CoMetaPNS). For reference, we further included a stationary MetaPNS model trained on all 12 tasks during meta-training. To assess the effect of meta-model pre-training on continual learning sample efficiency, we evaluated the continual learning methods under varying levels of data availability per new task, from 3 to 10 unique simulations. The stationary MetaPNS is trained with full data access, averaging 13.7 ± 2.5 samples per 26

Figure 7: Performance metrics evaluating the effectiveness of using a pre-trained metamodel as initialization for continual learning on held-out tasks when subject to varying levels of per-task sample availability.

task. Figure 7 summarizes results across all metrics, with a focus on RP for the continual learning methods. In the low-data regime, CoMetaPNS-pretrained demonstrated a marked improvement over both CoMetaPNS from-scratch and pure meta-generalization, even with just a few samples per task. As data availability increases, both CoMetaPNS-pretrained and CoMetaPNS tend to converge toward the performance of the full stationary model. More specifically, under limited data, CoMetaPNS from-scratch consistently underperformed against the CoMetaPNS-pretrained counterpart in SCC and MSE, indicating challenges in learning general propagation dynamics from sparse observations. The SCC gap narrows with increasing data. Curiously, for DC, a notable inversion appears: CoMetaPNS fromscratch surpasses even the full stationary version. We attribute this to an increased capacity of the meta-model to specialize towards the tissue-specific properties of these tasks, though at the cost of reduced propagation fidelity. 4.5. Additional Benefits of Task-Relational CoMetaPNS 4.5.1. Novel Subject Identification We further analyzed the benefits of the proposed GMM-based task-relational meta-learning in Fig. 8, where the average log-likelihood of each task’s test samples, relative to their nearest cluster mean, shows that the model effectively distinguishes between known and unknown data sources over time. Prior to a new data source’s integration, its samples exhibit consistently low 27

Figure 8: GMM log-likelihoods of all tasks’ samples over the sequence, becoming known to the model when its index is in the sequence.

likelihoods across all clusters, while once incorporated, the model reliably re-identifies and clusters future samples from the same source. To quantify this benefit of task identification strategy compared to the Task-Aware Meta-Learning strategy, we compared the computational performance of the fast-adaptation meta-model for a given identified known task versus always updating the weights of the meta-model. For the same number of iterations, data, and training structure - with the only difference being whether backpropagation was applied - we see an average speed improvement of 37.7% for processing all iterations. This provides a significant improvement in computational efficiency in the event of recurring, or similar, data distributions in the continual stream where we should not require retraining on known data. 4.5.2. Clustering Performance Fig. 9 illustrates the trajectory of context embeddings from the taskrelational reservoir across the task sequence, along with GMM cluster means, visualized in 2D via t-SNE. To align varying-sized embeddings resulting from different heart mesh geometries, shorter embeddings were padded with random noise to match the largest. While this introduces mesh-specific projection bias, clear intra-mesh clustering by material properties is still observed. 28

Figure 9: 2D t-SNE visualization of the task-relational reservoir over the task sequence. Black stars represent the GMM’s cluster means. X-marked samples are known. O-marked samples are unknown.

It can also be seen that there is clear separation between unknown and known tasks within the embedding space, supporting the use of log-likelihood as a means of novel subject identification. Given the set of identified known samples, for each underlying task ID, we compute the average purity of that task’s meta-embeddings under the GMM clustering. For each task, we identify all clusters containing at least one of its samples and calculate the purity of each, averaging these values to determine the task’s clustering purity. Fig. 10 shows these purities across the task sequence. While the GMM clustering achieves generally stable separation of task embeddings, it is not perfect, with an average final-task purity of 0.74 ± 0.18. Importantly, this analysis is based only on samples deemed "known" via log-likelihood thresholds, preserving the method’s ability to identify novel data sources. However, in some cases, the clustering fails to distinguish subtle differences in tissue properties on the same mesh structure. This can be seen in Fig. 9 where some of the cluster means across the sequence 29

Figure 10: Average cluster purity of each task’s samples over the task sequence. A purity of 1.0 indicates that all samples from a task are assigned to clusters containing only samples from that same task.

represent more than one task. 4.5.3. Effect of Rare Tasks As discussed previously, a key limitation of the Task-Aware Meta-Learning strategy is its assumption that each task corresponds to a distinct subject, which results in a fragmented and imbalanced reservoir as the number of task boundaries increases. To evaluate the impact of this on under-represented tasks, we conducted an ablation experiment with 25 subject-switching boundaries, where one task appeared only at the beginning and three other tasks cycled through the remaining boundaries. The reservoir was restricted to 200 samples. As shown in Fig. 11, the SCC performance of each subject over the sequence reveals that the Task-Aware strategy retained only about 5 samples from the first task by the end, leading to severe performance degradation. In contrast, the Task-Relational strategy maintained a more balanced allocation of approximately 50 samples per task, preserving performance even for the initially rare task.

30

Figure 11: Task-Relational vs. Task-Aware strategies on imbalanced tasks where the first task appears once while others cycle continually. Performance shown on the SCC metric.

4.5.4. Computational Requirements We evaluate the computational requirements of the GMM during training. This strategy introduces additional computational overhead only at task boundaries and does not affect adaptation efficiency at test time. GMM fitting is performed using the scikit-learn library on the CPU, requiring data transfer of meta-embeddings from the GPU, which introduces a minor slowdown. Compared to a meta-model implementation without GMM fitting, this addition resulted in only a 2.5% increase in training time, a modest cost given the benefits of task-relational modeling. Together, these results highlight the embedding space’s utility for robust subject-level identification at a limited computational cost. 5. Real-Data Experiments 5.1. Data, Baselines, & Experimental Setup Finally, we evaluated CoMetaPNS using in-vivo recordings from an animal model experiment (Bergquist et al., 2021), sourced from the Computational Electrocardiology Group. Cardiac activation sequences were generated via bipolar stimulation using intramural plunge needles at four distinct sites: left ventricular (LV) base, LV apex, LV freewall, and LV septum. Recordings were collected from a 247-electrode epicardial sock array with a spacing of 6.5 ± 1.3 mm. Geometric surfaces were reconstructed per recording, and five stimulations were performed at each location, yielding 20 total samples. We evaluated CoMetaPNS in two settings: 1) direct generalization of CoMetaPNS after it was continually trained on the synthetic data as described in Section 4, which we refer to as CoMetaPNS-MG (MG for metageneralization); and 2) pre-training of MetaPNS on the stationary multi-task 31

synthetic data as described in Section 4, and then continually fine-tuned on the real data, which we refer to CoMetaPNS-PT (PT for pre-trained). For comparison, we included baselines including the direct generalization of the stationary MetaPNS model and the continual, non-meta PNS. 5.2. Meta-Generalization of CoMetaPNS In a challenging meta-generalization setting, we tested the performance of the direct application of CoMetaPNS, MetaPNS, and PNS to real data when trained on synthetic data. For generalization testing, we selected one recording from a stimulation site as the target sample and used one recording from each of the remaining three locations as the context set. Cross-validation was conducted across all valid permutations, ensuring no overlap between target and context samples from the same stimulation site. Model performance on target samples was evaluated using mean squared error (MSE) and spatial correlation coefficient (SCC); Dice coefficient was omitted, as no scar tissue was present in these recordings. CoMetaPNS achieved an MSE of 8.9±0.9e-4 and SCC of 0.22±0.10, while MetaPNS recorded 8.3±0.04e-4 and 0.45±0.03, and PNS achieved 9.0±0.8e4 and 0.17±0.06. Figure 12 illustrates the mesh predictions from all three models. Although there is a notable performance drop compared to synthetic settings, CoMetaPNS continues to generalize better than the non-meta continual PNS. However, its gap to the stationary MetaPNS mirrors the generalization limitations observed in the synthetic experiments (Table 1). 5.3. CoMetaPNS with Pre-Training To assess whether pre-training accelerates the incorporation of novel information in real-world scenarios, we conducted additional experiments using the synthetic pre-trained model from Section 4.4 and adapted it to the in-vivo dataset through limited continual fine-tuning. For each experimental trial, one pacing site was designated as the held-out test set, while the remaining three sites constituted the training data. These training sites were partitioned into disjoint context-query sets per-batch and used for meta-learning over a fixed training window. Following this adaptation phase, we extracted a meta-embedding from the training pacing sites and evaluated the model’s performance on the held-out site, thereby measuring its capacity to acquire and transfer knowledge to unseen locations. To ensure statistical robustness, we performed leave-one-out cross-validation across all pacing sites and repeated each configuration across 5 random seeds. 32

Figure 12: Visual examples of reconstructed electrical activity from the proposed CoMetaPNS when generalizing from synthetic training data to in-vivo animal model experiment recordings in addition to continually fine-tuned results. Both epicardial and endocardial surfaces are presented. PT: Pre-Trained, MG: Meta-Generalization.

Compared to CoMetaPNS-MG when directly applied to the real data, CoMetaPNS-PT after continually fine-tuning on the real data achieved an MSE of 6.3±1.7e-4 and SCC of 0.56±0.08. These performance improvements are consistent with our previous findings when leveraging limited training data for rapid adaptation. Figure 12 illustrates the spatial predictions on the cardiac mesh, before and after fine-tuning. These results demonstrate that meta-models pre-trained exclusively on synthetic cardiac data can effectively bridge the domain gap to real cardiac recordings with minimal adaptation, despite the domain gap between computational simulations and in-vivo animal recordings. 6. Conclusion In this work, we introduced CoMetaPNS, a continual meta-learning framework designed to learn personalized neural surrogates from non-stationary streams of cardiac data. By combining a feed-forward meta-model with continual experience replay and Bayesian clustering for subject identification, CoMetaPNS enables efficient and scalable personalization without the need for retraining. Empirical results show that maintaining predictive performance across dynamically evolving subject distributions requires the joint use of both meta-learning and continual learning components. This frame33

work provides a practical approach for real-world, personalized cardiac simulation, and future work will explore its extension to diverse observational data types and more complex clinical applications. This study was conducted extensively on synthetic data with a feasibility study on real animal-model data. This constraint is primarily due to the authors’ lack of access to in-vivo recordings from experimental or clinical subjects undergoing cardiac mapping and stimulation procedures. The feasibility study of the ability of CoMetaPNS-PT to quickly adapt to real data using samples from only three pacing sites however provides a glimpse of the significant potential for CoMetaPNS-PT to be used with limited real data. Future work will focus on this aspect. Various modeling choices exist for representing cardiac electrical activity in virtual heart simulations, including heart surface potentials (Malik et al., 2018; Karoui et al., 2019; Horvath et al., 2019) and transmembrane voltages defined over the volumetric mesh of the heart (Wang et al., 2009; Dhamala et al., 2019; Ghimire et al., 2018, 2019a,b). This study focuses on the former, as it represents the most widely adopted formulation and is used in commercial systems. Extending CoMetaPNS to volumetric representations presents additional challenges, particularly in constructing suitable hierarchical graph structures for spatial decoding. Moreover, the significantly larger graph size associated with volumetric models would further increase the computational cost of training. A complete cardiac signal consists of two primary stages: depolarization and repolarization. However, accurately determining the timing of local electrical activation and recovery from transmembrane voltages or extracellular potentials remains an open challenge. In this study, we restricted our focus to the depolarization stage due to the current limitations of neural networks in capturing more complex physiological dynamics. Future work will explore the use of more advanced neural architectures to better model the full sequence of cardiac electrical activity. 7. Declaration of generative AI and AI-assisted technologies in the manuscript preparation process. During the preparation of this work, the author(s) used Claude Sonnet 4.5 in order to perform writing style and grammar tweaks. No generative tools were used for research purposes, literature review, or content generation. After using this tool/service, the author(s) reviewed and edited the content 34

as needed and take(s) full responsibility for the content of the published article. Acknowledgments We thank the Computational Electrocardiology Group for their support in providing access to and resources for the in-vivo animal heart data used within Section 5. This study was supported by National Institutes of Health grant 1R01HL145590, 2R01HL145590 and R01NR018301, and National Science Foundation grant number OAC-2212548. References Aliev, R.R., Panfilov, A.V., 1996. A simple two-variable model of cardiac excitation. Chaos, Solitons & Fractals 7, 293–301. Arevalo, H.J., Vadakkumpadan, F., Guallar, E., Jebb, A., Malamas, P., Wu, K.C., Trayanova, N.A., 2016. Arrhythmia risk stratification of patients after myocardial infarction using personalized heart models. Nature communications 7, 1–8. Bergquist, J., Good, W., Zenger, B., Tate, J., Rupp, L., MacLeod, R., 2021. The electrocardiographic forward problem: A benchmark study. Comput Biol Med 134, 104476. Bhagirath, P., Strocchi, M., Bishop, M.J., Boyle, P.M., Plank, G., 2024. From bits to bedside: entering the age of digital twins in cardiac electrophysiology. Europace 26, euae295. Brebbia, C.A., Telles, J.C.F., Wrobel, L.C., 2012. Boundary element techniques: theory and applications in engineering. Springer Science & Business Media. Caccia, M., Rodriguez, P., Ostapenko, O., Normandin, F., Lin, M., PageCaccia, L., Laradji, I.H., Rish, I., Lacoste, A., Vázquez, D., et al., 2020. Online fast adaptation and knowledge accumulation (osaka): a new approach to continual learning. Advances in Neural Information Processing Systems 33, 16532–16545.

35

Cacciola, F., 2007. Triangulated surface mesh simplification, in: Board, C.E. (Ed.), CGAL User and Reference Manual. 3.3 ed. URL: http://www.cgal.org/Manual/3.3/doc_html/cgal_manual/packages. html#Pkg:SurfaceMeshSimplification. Cantwell, C.D., Mohamied, Y., Tzortzis, K.N., Garasto, S., Houston, C., Chowdhury, R.A., Ng, F.S., Bharath, A.A., Peters, N.S., 2019. Rethinking multiscale cardiac electrophysiology with machine learning and predictive modelling. Computers in biology and medicine 104, 339–351. Chung, J., Gulcehre, C., Cho, K., Bengio, Y., 2014. Empirical evaluation of gated recurrent neural networks on sequence modeling. arXiv preprint arXiv:1412.3555 . Cluitmans, M.J., Plank, G., Heijman, J., 2024. Digital twins for cardiac electrophysiology: state of the art and future challenges. Herzschrittmachertherapie+ Elektrophysiologie 35, 118–123. Coveney, S., Corrado, C., Oakley, J.E., Wilkinson, R.D., Niederer, S.A., Clayton, R.H., 2021. Bayesian calibration of electrophysiology models using restitution curve emulators. Frontiers in Physiology 12, 1120. Dawoud, F., 2009. Noninvasive Imaging of Epicardial Potentials for Clinical Electrophysiology. Ph.D. thesis. Dalhousie University. Dhamala, J., Arevalo, H.J., Sapp, J., Horácek, B.M., Wu, K.C., Trayanova, N.A., Wang, L., 2018a. Quantifying the uncertainty in model parameters using gaussian process-based markov chain monte carlo in cardiac electrophysiology. Medical image analysis 48, 43–57. Dhamala, J., Bajracharya, P., Arevalo, H.J., Sapp, J.L., Horácek, B.M., Wu, K.C., Trayanova, N.A., Wang, L., 2020. Embedding high-dimensional bayesian optimization via generative modeling: parameter personalization of cardiac electrophysiological models. Medical image analysis 62, 101670. Dhamala, J., Ghimire, S., Sapp, J.L., Horáček, B.M., Wang, L., 2018b. Highdimensional bayesian optimization of personalized cardiac model parameters via an embedded generative model, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 499–507. 36

Dhamala, J., Ghimire, S., Sapp, J.L., Horáček, B.M., Wang, L., 2019. Bayesian optimization on large graphs via a graph convolutional generative model: Application in cardiac model personalization, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 458–467. Fey, M., Eric Lenssen, J., Weichert, F., Müller, H., 2018. Splinecnn: Fast geometric deep learning with continuous b-spline kernels, in: The IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pp. 869–877. Finn, C., Abbeel, P., Levine, S., 2017a. Model-agnostic meta-learning for fast adaptation of deep networks, in: International conference on machine learning, PMLR. pp. 1126–1135. Finn, C., Abbeel, P., Levine, S., 2017b. Model-agnostic meta-learning for fast adaptation of deep networks, in: International conference on machine learning, PMLR. pp. 1126–1135. Fresca, S., Manzoni, A., Dedè, L., Quarteroni, A., 2021. Pod-enhanced deep learning-based reduced order models for the real-time simulation of cardiac electrophysiology in the left atrium. Frontiers in physiology 12, 1431. Ghimire, S., Dhamala, J., Gyawali, P.K., Sapp, J.L., Horacek, M., Wang, L., 2018. Generative modeling and inverse imaging of cardiac transmembrane potential, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 508–516. Ghimire, S., Gyawali, P.K., Dhamala, J., Sapp, J.L., Horacek, M., Wang, L., 2019a. Improving generalization of deep networks for inverse reconstruction of image sequences, in: International Conference on Information Processing in Medical Imaging, Springer. pp. 153–166. Ghimire, S., Sapp, J.L., Horáček, B.M., Wang, L., 2019b. Noninvasive reconstruction of transmural transmembrane potential with simultaneous estimation of prior model error. IEEE transactions on medical imaging 38, 2582–2595. Giffard-Roisin, S., Jackson, T., Fovargue, L., Lee, J., Delingette, H., Razavi, R., Ayache, N., Sermesant, M., 2016. Noninvasive personalization of a car37

diac electrophysiology model from body surface potential mapping. IEEE Transactions on Biomedical Engineering 64, 2206–2218. Gomez, R.A., Stöcker, J., Cansız, B., Kaliske, M., 2025. Simulation of parametrized cardiac electrophysiology in three dimensions using physicsinformed neural networks. arXiv preprint arXiv:2506.15405 . Harrison, J., Sharma, A., Finn, C., Pavone, M., 2020. Continuous metalearning without tasks. Advances in neural information processing systems 33, 17571–17581. He, X., Sygnowski, J., Galashov, A., Rusu, A.A., Teh, Y.W., Pascanu, R., 2019. Task agnostic continual learning via meta learning. arXiv preprint arXiv:1906.05201 . Herrero Martin, C., Oved, A., Chowdhury, R.A., Ullmann, E., Peters, N.S., Bharath, A.A., Varela, M., 2022. Ep-pinns: Cardiac electrophysiology characterisation using physics-informed neural networks. Frontiers in Cardiovascular Medicine 8, 768419. Horvath, J., Shien, L., Peng, T., Malik, A., Trew, M., Bear, L., 2019. Deep learning neural nets for detecting heart activity. arXiv preprint arXiv:1901.09831 . Jiang, X., Ghimire, S., Dhamala, J., Li, Z., Gyawali, P.K., Wang, L., 2020. Learning geometry-dependent and physics-based inverse image reconstruction, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 487–496. Jiang, X., Li, Z., Missel, R., Zaman, M.S., Zenger, B., Good, W.W., MacLeod, R.S., Sapp, J.L., Wang, L., 2022. Few-shot generation of personalized neural surrogates for cardiac simulation via bayesian meta-learning, in: International Conference on Medical Image Computing and ComputerAssisted Intervention, Springer. pp. 46–56. Jiang, X., Missel, R., Toloubidokhti, M., Li, Z., Gharbia, O., Sapp, J.L., Wang, L., 2021. Label-free physics-informed image sequence reconstruction with disentangled spatial-temporal modeling, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 361–371. 38

Joseph, J., Gu, A., 2021. La-maml: Look-ahead meta learning for continual learning, ml reproducibility challenge 2020. arXiv preprint arXiv:2102.05824 . Karoui, A., Bendahmane, M., Zemzemi, N., 2019. A spatial adaptation of the time delay neural network for solving ecgi inverse problem, in: International Conference on Functional Imaging and Modeling of the Heart, Springer. pp. 94–102. Kashtanova, V., Ayed, I., Cedilnik, N., Gallinari, P., Sermesant, M., 2021. Ep-net 2.0: out-of-domain generalisation for deep learning models of cardiac electrophysiology, in: International Conference on Functional Imaging and Modeling of the Heart, Springer. pp. 482–492. Kasim, M.F., Watson-Parris, D., Deaconu, L., Oliver, S., Hatfield, H., Froula, D.H., Gregori, G., Jarvis, M., Khatiwala, S., Korenaga, J., ToppMugglestone, J., Viezzer, E., Vinko, S.M., 2021. Building high accuracy emulators for scientific simulations with deep neural architecture search. Machine Learning: Science and Technology , in press. Kingma, D.P., Ba, J., 2014. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 . Kingma, D.P., Welling, M., 2013. Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114 . Kirkpatrick, J., Pascanu, R., Rabinowitz, N., Veness, J., Desjardins, G., Rusu, A.A., Milan, K., Quan, J., Ramalho, T., Grabska-Barwinska, A., et al., 2017. Overcoming catastrophic forgetting in neural networks. Proceedings of the national academy of sciences 114, 3521–3526. Kolk, M.Z., Ruipérez-Campillo, S., Alvarez-Florez, L., Deb, B., Bekkers, E.J., Allaart, C.P., Van Der Lingen, A.L.C., Clopton, P., Išgum, I., Wilde, A.A., et al., 2024. Dynamic prediction of malignant ventricular arrhythmias using neural networks in patients with an implantable cardioverterdefibrillator. EBioMedicine 99. Liu, G.R., 2003. Mesh free methods: moving beyond the finite element method. CRC Press, Boca Raton, Fla.

39

Malik, A., Peng, T., Trew, M.L., 2018. A machine learning approach to reconstruction of heart surface potentials from body surface potentials, in: 2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), IEEE. pp. 4828–4831. Martinez, E., Moscoloni, B., Salvador, M., Kong, F., Peirlinck, M., Lesley Marsden, A., 2025. Full-field surrogate modeling of cardiac function encoding geometric variability . Miller, R., Kerfoot, E., Mauger, C., Ismail, T.F., Young, A.A., Nordsletten, D.A., 2021. An implementation of patient-specific biventricular mechanics simulations with a deep learning and computational pipeline. Frontiers in physiology 12, 1398. Mishra, N., Rohaninejad, M., Chen, X., Abbeel, P., 2017. A simple neural attentive meta-learner. arXiv preprint arXiv:1707.03141 . Neumann, D., Mansi, T., 2020. Machine learning methods for robust parameter estimation, in: Artificial Intelligence for Computational Modeling of the Heart. Elsevier, pp. 161–181. Niederer, S., Aboelkassem, Y., Cantwell, C.D., Corrado, C., Coveney, S., Cherry, E.M., Delhaas, T., Fenton, F.H., Panfilov, A., Pathmanathan, P., et al., 2020. Creation and application of virtual patient cohorts of heart models. Philosophical Transactions of the Royal Society A 378, 20190558. Otsu, N., 1979. A threshold selection method from gray-level histograms. IEEE transactions on systems, man, and cybernetics 9, 62–66. Prakosa, A., Arevalo, H.J., Deng, D., Boyle, P.M., Nikolov, P.P., Ashikaga, H., Blauer, J.J., Ghafoori, E., Park, C.J., Blake, R.C., et al., 2018. Personalized virtual-heart technology for guiding the ablation of infarct-related ventricular tachycardia. Nature biomedical engineering 2, 732–740. Ravi, S., Larochelle, H., 2016. Optimization as a model for few-shot learning, in: International conference on learning representations. Riemer, M., Cases, I., Ajemian, R., Liu, M., Rish, I., Tu, Y., Tesauro, G., 2018. Learning to learn without forgetting by maximizing transfer and minimizing interference. arXiv preprint arXiv:1810.11910 . 40

Salvador, M., Strocchi, M., Regazzoni, F., Augustin, C.M., Dede’, L., Niederer, S.A., Quarteroni, A., 2024. Whole-heart electromechanical simulations using latent neural ordinary differential equations. NPJ Digital Medicine 7, 90. Sermesant, M., Chabiniok, R., Chinchapatnam, P., Mansi, T., Billet, F., Moireau, P., Peyrat, J.M., Wong, K., Relan, J., Rhode, K., et al., 2012. Patient-specific electromechanical models of the heart for the prediction of pacing acute effects in crt: a preliminary clinical validation. Medical image analysis 16, 201–215. Trayanova, N.A., Lyon, A., Shade, J., Heijman, J., 2024. Computational modeling of cardiac electrophysiology and arrhythmogenesis: toward clinical translation. Physiological reviews 104, 1265–1333. Wang, L., Zhang, H., Wong, K.C., Liu, H., Shi, P., 2009. Physiologicalmodel-constrained noninvasive reconstruction of volumetric myocardial transmembrane potentials. IEEE Transactions on Biomedical Engineering 57, 296–315. Wong, K.C., Sermesant, M., Rhode, K., Ginks, M., Rinaldi, C.A., Razavi, R., Delingette, H., Ayache, N., 2015. Velocity-based cardiac contractility personalization from images using derivative-free optimization. Journal of the mechanical behavior of biomedical materials 43, 35–52. Zettinig, O., Mansi, T., Georgescu, B., Kayvanpour, E., Sedaghat-Hamedani, F., Amr, A., Haas, J., Steen, H., Meder, B., Katus, H., et al., 2013. Fast data-driven calibration of a cardiac electrophysiology model from images and ecg, in: International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. pp. 1–8.

41