Why Cardiac Digital Twins Can’t Be Solved From One ECG

Analysis by the aitrendblend editorial team · Pillar 1, AI for Medical Imaging and Healthcare · 15 minute read

Picture a cardiologist staring at a 12 lead ECG trace from a patient being considered for cardiac resynchronization therapy, the kind of device that paces both sides of the heart back into rhythm. The tracing tells her something is off. What it does not tell her, not directly, is exactly how electricity is spreading through this particular patient’s ventricles, which is precisely the detail that would let a surgeon target the pacing leads with confidence instead of an educated guess. A team from the University of Graz, Universita della Svizzera italiana, and the University of Trento just built a tool that tries to answer that question from the ECG alone. What they found along the way is almost as important as the tool itself.

The part of the heart you cannot see on a standard scan

Every heartbeat starts with an electrical signal racing through a specialized wiring system called the His Purkinje system, or HPS. It begins near the atrioventricular node, runs down the His Bundle, splits into left and right bundle branches, and fans out into a fast conducting Purkinje network that hugs the inner surface of both ventricles. Where that network actually touches ventricular muscle, at points called Purkinje muscle junctions or PMJs, the electrical wave jumps from the wiring into the muscle itself and triggers the contraction that pumps blood out to the body. Get this wiring diagram right for a specific patient and you can, in principle, simulate exactly how their heart will respond to a pacemaker lead, a drug, or a surgical ablation before you ever touch the patient.

The frustrating part is that this wiring is essentially invisible. It cannot be resolved with standard MRI or CT. Even invasive catheter mapping, threading electrodes directly into the heart chambers, only samples it with limited accuracy (Palamara et al., 2014). So researchers building cardiac digital twins, which the authors define as one to one computational replicas of a specific patient’s heart, are stuck trying to infer this hidden structure from whatever they can measure safely from outside the body, which usually means the 12 lead ECG (Niederer et al., 2021 and Corral Acero et al., 2020).

Why most digital twins never actually get personalized

Here is an uncomfortable admission buried in the introduction of this paper. Building a patient specific heart shape from medical images, what the field calls anatomical twinning, has become fairly automated (Crozier et al., 2016). Functional twinning, tuning the electrical parameters so the model’s simulated ECG actually matches the patient’s real ECG, has lagged far behind. Most published digital twin studies quietly sidestep this step. They plug in the same population average electrical parameters for every patient and never do a like for like comparison against the patient’s real recorded signal, because doing so would expose how far off the model actually is (Sung et al., 2022 and Bishop and Plank, 2025).

That gap matters enormously if you take the word twin seriously. A model that was never checked against the patient’s own ECG is not really personalized, it is a generic heart shape wearing the patient’s anatomy. Closing that gap is what this paper is actually about, and the tool they built to do it is called Geodesic-BP, first introduced in earlier work by the same group (Grandits et al., 2024) and put through a much more thorough stress test here.

Teaching backpropagation to walk geodesic paths through heart muscle

The forward model at the center of all this is the eikonal equation, a standard way of simulating how an electrical wavefront spreads through heart tissue without having to solve the full, much more expensive cellular scale reaction diffusion equations.

Once you know the activation time at every point in the ventricles, you can reconstruct the transmembrane voltage as a smooth wavefront passing through that point, then use the lead field method to translate that voltage pattern into what each of the 12 ECG leads would actually record.

None of this is new on its own. Fast eikonal solvers have existed for years (Neic et al., 2017 and Pezzuto et al., 2017). What makes Geodesic-BP different is treating the whole pipeline, from PMJ positions and timings all the way through to the simulated ECG, as one differentiable computation, and then using standard machine learning backpropagation to compute the gradient of the mismatch with respect to every PMJ. That gradient tells the optimizer exactly which direction to nudge each junction to better match the recorded ECG. The name comes from the fact that backpropagating through the fixed point eikonal solver effectively returns discrete geodesic paths from each PMJ, the shortest travel time routes through the anisotropic tissue.

The optimization itself runs for 400 iterations of the Adam optimizer (Kingma and Ba, 2017) with a learning rate of 0.75, and the authors report that a single iteration now takes about 6 seconds after code optimization, down from roughly 16 seconds in the original implementation. On a single high end NVidia RTX 4090 GPU, fitting 300 PMJs to a full 12 lead ECG takes under 30 minutes. The same job on a CPU stretches past 6 hours, which is the difference between a method a cardiac electrophysiology lab could actually run overnight and one that stays a research curiosity.

Building anatomy into the search space

Letting the optimizer place PMJs absolutely anywhere in the ventricular muscle sounds like the more flexible choice, and in one sense it is. The paper shows it is also a recipe for anatomically nonsensical answers. So the team also tested a restricted version, where PMJs are only allowed to sit inside a subendocardial band, a thin shell hugging the inner heart surface, matching where the real Purkinje network actually lives.

That 2.5 millimeter cutoff was not reverse engineered from the specific ground truth model used for testing. It came from general population data on ventricular wall thickness (Sjogren, 1971 and Prakash, 1978) and from histological evidence that Purkinje fibers rarely penetrate more than about 30 percent of the wall depth (Vigmond and Stuyvers, 2016). The authors are explicit that they deliberately avoided peeking at the exact geometry of their own test heart when setting this constraint, specifically so the result would say something about physiological priors in general rather than about overfitting to one test case.

What the ground truth heart actually was

To check whether any of this actually works, you need a case where you know the real answer, which is nearly impossible with a living patient. The authors instead built a synthetic ground truth from MRI scans of a single healthy male subject, comprising a torso scan and a whole heart cardiac scan, with real ECG recordings from ten electrodes taken during the same imaging session. A convolutional neural network handled automatic cardiac segmentation from the MRI (Payer et al., 2018), and a finite element mesh was generated with roughly 1253 micrometers of resolution in the ventricles.

Onto that anatomy the team grafted a highly detailed, topologically realistic Purkinje network containing 907 real PMJs, discretized at a spatial resolution of 482 micrometers, calibrated in earlier work to reproduce this specific subject’s ECG with high fidelity (Gillette et al., 2022). This detailed model became the answer key. Geodesic-BP never got to see its topology or its true PMJ count. It only ever saw the resulting simulated ECG and had to work backward from that alone.

The result nobody wanted to find, and could not ignore

Run Geodesic-BP once on this ECG and the fit looks spectacular. The maximum error between the optimized and ground truth ECG was under 0.028 millivolts, with an average error around 0.0207 millivolts and a Pearson correlation above 0.994 across all runs. Convergence to within 0.1 millivolts typically happened in 100 iterations or fewer. On paper, this is about as clean a fit as you could hope for.

Then the authors did something most papers in this space skip entirely. They ran the same fitting procedure 20 separate times, starting each run from a different random guess of 300 PMJ positions and timings, and compared the resulting internal activation maps to each other, not just to the ECG. The ECGs from all 20 runs still looked nearly identical to the target. The underlying activation maps, the actual wiring the model inferred, did not. The two most divergent runs differed by an average of 18.32 milliseconds across the whole ventricle, with a peak local difference of 23.12 milliseconds. In one run, a spot on the mid anterior wall of the left ventricle was among the very first places to activate. In another run fit to the exact same ECG, that same spot was among the very last.

A tight fit to the ECG is not proof of a correct answer underneath it. Two very different hearts can sound identical from the outside. Editorial synthesis of the study’s central identifiability finding

Body surface potential maps told a similar story. Using a denser 128 electrode torso vest, the two most different reconstructed potential maps differed by up to 0.11 millivolts at their peak, and at one unseen point on the sternum the correlation between two competing solutions and the true potential dropped to just 0.69 and 0.64. The overall dipole shaped pattern of the ECG was preserved everywhere. The fine detail, especially near the chest wall directly over the heart, was not.

Physiology closes the gap that more electrodes cannot

Given that unrestricted fitting produces multiple, meaningfully different answers, the obvious next question is what actually helps. The paper tests two candidate fixes against each other, physiological constraints on where PMJs can sit, and simply recording from more electrodes.

Read down each column and the pattern jumps out. Going from a 12 lead ECG to a 128 electrode vest, a huge jump in the amount of data collected, only nudged the unrestricted error from 20.08 down to 18.95 milliseconds. Adding the anatomical constraint, with no extra hardware at all, dropped the 12 lead result all the way from 20.08 to 14.63 milliseconds, a far larger improvement from a far cheaper intervention. Averaged across everything, the physiologically restricted models settled to an activation timing spread of 8.38 to 12.13 milliseconds, compared with 12.13 to over 20 milliseconds unrestricted, and the authors note this spread is already narrower than the accuracy typically claimed for the best invasive electro anatomical mapping catheters used in real procedures today.

How many starting guesses does the optimizer actually need

One practical question any lab trying to reuse this method will ask immediately is how many candidate PMJs to start with. Too few and the model simply cannot represent the activation pattern. The authors tested counts from 1 up to 5000. Below about 20 initial PMJs, the fit failed outright, with ECG correlation under 0.8 and activation timing errors above 30 milliseconds. Between 100 and 500 PMJs, the ECG fit essentially saturated, with correlations of 0.998 and near identical reconstructions. Push past 500 and the ECG fit improved only marginally while the optimization became noticeably more sensitive to the random starting positions.

There is also a quiet efficiency story hiding in these numbers. With an initial pool of 300 PMJs, only about 23 percent of them, roughly 70 junctions, ended up doing 95 percent of the work of activating the ventricle by the end of optimization. Push the initial pool to 5000 and that active fraction fell to about 11 percent, close to 545 sites. The rest of the disabled candidates were not wasted computation exactly, they acted more like an oversized search net that the optimizer was free to prune, which is a very different failure mode than classical sampling methods that choke as parameter counts grow.

Clinical translation gap

It is worth being direct about the distance between this result and anything a hospital would use tomorrow. Everything in this paper, the fitting accuracy, the identifiability findings, and the effect of physiological constraints, was measured against a single synthetic ground truth heart, built from one healthy male volunteer’s MRI and modeled with a separate, independently developed simulator called CARPentry (Vigmond et al., 2008). That is a controlled and genuinely clever way to get an answer key that would otherwise be impossible to obtain noninvasively, but it means the entire quantitative story rests on one anatomy, one healthy rhythm, and one set of assumed conduction velocities.

The paper has not yet been tested against a patient with an actual conduction disorder, a bundle branch block, or a diseased or scarred ventricle, all situations where the anatomy and electrical behavior look nothing like the healthy test case used here. It also stops at the QRS complex, the depolarization phase of the heartbeat, and explicitly leaves out the repolarization phase and the T wave, which carries its own clinically important information about arrhythmia risk. And the method as described requires a high end consumer GPU to hit its 30 minute runtime, a specific hardware dependency that any clinical deployment would need to account for. None of this diminishes the engineering achievement. It does mean the honest way to describe this work is as a strong proof of concept for the identifiability problem in cardiac twinning, not as a validated clinical tool.

Honest limitations, with the paper’s own numbers

Where this could go next

The most interesting open door the authors leave is a Bayesian reframing of the whole problem. Right now Geodesic-BP returns a point estimate, plus an empirical spread from running it many times with different starting guesses. A formal Bayesian treatment, with a genuine prior distribution over plausible Purkinje network shapes built from rule based or generative anatomical models, could turn that empirical spread into a properly calibrated uncertainty estimate a clinician could actually reason about. The authors also flag extending the fitting objective to the repolarization phase and the T wave, optimizing the conduction velocity tensor itself rather than holding it fixed, and building denser anatomical models that keep structures like the moderator band and papillary muscles as actual tissue rather than folding them into the blood pool.

Reproducing the core optimization loop

The full Geodesic-BP pipeline depends on a finite element mesh solver, GPU accelerated eikonal fast marching, and precomputed lead fields from a torso model, which is far more infrastructure than fits cleanly into a single teaching example. The implementation below distills the mathematical core the paper actually optimizes, representing the ventricular domain as a point cloud, approximating the eikonal travel time with a differentiable soft minimum over pairwise geodesic style distances, and fitting PMJ positions and timings to a target ECG using PyTorch and the Adam optimizer, exactly as described in Section 2.3 of the paper. It is a simplified stand in for the mesh based solver, built so the fitting loop, loss function, and physiological constraint logic can be tested and understood without a full cardiac mesh pipeline.

# geodesic_bp_lite.py
# A simplified, mesh free reproduction of the Geodesic-BP fitting loop
# from Grandits, Gillette, Plank and Pezzuto, Medical Image Analysis 105 (2025) 103641.
# Approximates the eikonal activation time with a differentiable soft minimum
# over travel times from candidate PMJs, then fits PMJ position and timing
# to a target ECG using gradient descent, matching Eq. 10 of the paper.

import torch
import torch.nn as nn


class GeodesicBPLite(nn.Module):
    """Differentiable stand in for the eikonal forward model and lead field
    ECG computation. Ventricular tissue is a point cloud, PMJs are trainable
    positions and activation times, and the ECG is a linear projection of
    the resulting soft minimum activation map."""

    def __init__(self, mesh_points: torch.Tensor, lead_field: torch.Tensor,
                 num_pmjs=300, conduction_velocity=0.6, softness=8.0):
        super().__init__()
        self.mesh_points = mesh_points          # (n_nodes, 3), fixed ventricular geometry
        self.lead_field = lead_field            # (n_leads, n_nodes), fixed precomputed weights
        self.cv = conduction_velocity            # meters per second, treated as known
        self.softness = softness

        n_nodes = mesh_points.shape[0]
        init_idx = torch.randperm(n_nodes)[:num_pmjs]
        self.pmj_positions = nn.Parameter(mesh_points[init_idx].clone())
        self.pmj_times = nn.Parameter(torch.rand(num_pmjs) * 10.0)  # ms

    def project_to_feasible_region(self, feasible_center: torch.Tensor, feasible_radius: float):
        """Implements the restricted case from Section 2.5, clamping each PMJ
        to lie within feasible_radius of the endocardial reference points,
        and clamping activation times to be non negative as in Eq. 11."""
        with torch.no_grad():
            offset = self.pmj_positions - feasible_center
            dist = torch.norm(offset, dim=1, keepdim=True)
            too_far = (dist > feasible_radius).squeeze()
            if too_far.any():
                scale = feasible_radius / torch.clamp(dist[too_far], min=1e-6)
                self.pmj_positions[too_far] = feasible_center + offset[too_far] * scale
            self.pmj_times.clamp_(min=0.0)

    def activation_map(self) -> torch.Tensor:
        """Soft minimum approximation of the eikonal LAT, Eq. 1 and Eq. 5.
        A true fast marching solver enforces an exact minimum. Here a
        soft minimum keeps the whole computation differentiable end to end."""
        diffs = self.mesh_points.unsqueeze(1) - self.pmj_positions.unsqueeze(0)  # (n_nodes, n_pmj, 3)
        travel_dist = torch.norm(diffs, dim=-1)                              # (n_nodes, n_pmj)
        candidate_times = self.pmj_times.unsqueeze(0) + travel_dist / self.cv
        # differentiable soft minimum instead of a hard min, per node
        weights = torch.softmax(-self.softness * candidate_times, dim=1)
        lat = (weights * candidate_times).sum(dim=1)
        return lat  # (n_nodes,)

    def transmembrane_voltage(self, lat: torch.Tensor, t: float, v0=-85.0, v1=30.0, eps=1.0) -> torch.Tensor:
        """Bistable upstroke model, Eq. 6 of the paper."""
        return v0 + (v1 - v0) / 2.0 * (torch.tanh(2.0 * (t - lat) / eps) + 1.0)

    def forward(self, time_grid: torch.Tensor) -> torch.Tensor:
        """Returns the simulated ECG, shape (n_leads, n_timesteps), Eq. 8."""
        lat = self.activation_map()
        ecg = []
        for t in time_grid:
            vm = self.transmembrane_voltage(lat, t.item())
            ecg.append(self.lead_field @ vm)
        return torch.stack(ecg, dim=1)


def fit_geodesic_bp(model: GeodesicBPLite, target_ecg: torch.Tensor, time_grid: torch.Tensor,
                     feasible_center=None, feasible_radius=None,
                     iterations=400, lr=0.75):
    """Gradient based fitting loop matching Section 2.3. Uses Adam for
    400 iterations at learning rate 0.75, and projects PMJs back onto
    the feasible set after every step when constraints are supplied."""
    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
    loss_history = []

    for iteration in range(iterations):
        optimizer.zero_grad()
        simulated_ecg = model(time_grid)
        loss = torch.mean((simulated_ecg - target_ecg) ** 2)
        loss.backward()
        optimizer.step()

        if feasible_center is not None:
            model.project_to_feasible_region(feasible_center, feasible_radius)

        loss_history.append(loss.item())

    return loss_history


def dist_v(simulated: torch.Tensor, target: torch.Tensor, dt: float) -> float:
    """Root mean squared ECG mismatch, Eq. 17 of the paper."""
    n_leads, n_t = simulated.shape
    mse = torch.mean((simulated - target) ** 2)
    return torch.sqrt(mse).item()


def run_smoke_test():
    """Builds a dummy ventricle as a point cloud, a random ground truth
    PMJ set, and confirms Geodesic-BP Lite can recover an ECG match
    starting from a different random PMJ initialization."""
    torch.manual_seed(0)
    n_nodes, n_leads, n_pmj_gt = 800, 12, 40

    mesh_points = torch.randn(n_nodes, 3) * 40.0            # mm scale ventricle point cloud
    lead_field = torch.randn(n_leads, n_nodes) * 0.01       # stand in for a precomputed lead field
    time_grid = torch.linspace(0.0, 120.0, 60)             # ms, covers a QRS complex

    # ground truth model, held fixed, used only to generate a target ECG
    gt_model = GeodesicBPLite(mesh_points, lead_field, num_pmjs=n_pmj_gt)
    with torch.no_grad():
        target_ecg = gt_model(time_grid)

    # fitting model, randomly initialized, restricted to a subendocardial style shell
    fit_model = GeodesicBPLite(mesh_points, lead_field, num_pmjs=300)
    feasible_center = mesh_points.mean(dim=0)
    feasible_radius = 45.0  # mm, generous stand in for the 2.5 mm endocardial band concept

    loss_history = fit_geodesic_bp(
        fit_model, target_ecg, time_grid,
        feasible_center=feasible_center, feasible_radius=feasible_radius,
        iterations=400, lr=0.75,
    )

    with torch.no_grad():
        final_ecg = fit_model(time_grid)
        final_error = dist_v(final_ecg, target_ecg, dt=time_grid[1].item() - time_grid[0].item())

    print(f"Initial loss {loss_history[0]:.5f}, final loss {loss_history[-1]:.5f}")
    print(f"Final ECG RMS error {final_error:.5f}")

    assert loss_history[-1] < loss_history[0]
    assert final_ecg.shape == target_ecg.shape
    print("Smoke test passed.")


if __name__ == "__main__":
    run_smoke_test()

The bigger picture for digital twins in medicine

Step back from the millisecond by millisecond results and this paper is making a point that applies well beyond cardiology. Every field racing toward personalized digital twins, whether of a heart, a lung, or a tumor, faces the same quiet risk. A model can match your measured data beautifully and still be wrong about the hidden structure that produced that data, simply because more than one hidden structure can produce it. Fitting error alone cannot tell you which case you are in. You need either more informative measurements, a genuine physiological prior, or an honest accounting of uncertainty, and this paper is a rare example of a group actually measuring which of those levers matters most rather than assuming the answer.

The conceptual shift worth carrying forward is that identifiability deserves to be treated as a first class question, not an afterthought mentioned in a limitations paragraph. The authors did not just build a faster fitting algorithm. They used that speed to run the same fit hundreds of times and interrogate whether the fitting problem even has a stable answer, which is a genuinely different scientific move than reporting a single best fit and calling it a day.

The transferability angle is real but should not be overstated. The gradient based, differentiable forward model approach could plausibly extend to other biophysical inverse problems built on similarly structured PDEs, and the general lesson about anatomical priors beating denser sensors is intuitive enough to test in adjacent fields. But the specific numbers in this paper, the millisecond errors and the electrode counts, are tied tightly to cardiac electrophysiology and to this one synthetic test case, and porting the method elsewhere would require its own validation work from scratch.

What remains genuinely unresolved is the leap from this synthetic, single subject proof of concept to a real diseased heart with scar tissue, an accessory pathway, or an existing pacemaker, none of which resemble the clean healthy anatomy tested here. The authors are candid about this gap themselves, flagging repolarization modeling, joint estimation of conduction velocity, and Bayesian priors over Purkinje anatomy as the next steps rather than finished work.

Read plainly, this is a methods paper that earns its place by being honest about a hard problem rather than by claiming to have solved it. It hands the field a faster tool, a sharper way to quantify how uncertain a fitted cardiac twin actually is, and a clear demonstration that biology, not more electrodes, is the more efficient path to trustworthy answers. That combination, speed paired with epistemic honesty, is a rarer thing in computational medicine than it should be.

Frequently asked questions

Explore more from this pillar