Introduction
The human heart is an extraordinarily complex organ, and understanding its electrical behavior has long been one of medicine’s greatest challenges. For patients suffering from atrial fibrillation (AF) and other cardiac rhythm disorders, traditional treatment approaches rely heavily on trial-and-error methodologies and preclinical animal testing. However, a revolutionary breakthrough in cardiac imaging and computational modeling is changing everything: the creation of digital twin snapshots—virtual replicas of individual patient hearts that can predict how their unique anatomy and physiology will respond to therapeutic interventions before a single procedure is performed.
A groundbreaking study published in Medical Image Analysis describes a comprehensive automated workflow that represents a critical leap forward in creating patient-specific digital twins of human atria with unprecedented anatomical accuracy and computational efficiency. This advancement promises to transform how clinicians approach atrial fibrillation ablation, device development, and personalized cardiac therapy planning.
Understanding Digital Twins in Cardiac Medicine
What Are Cardiac Digital Twins?
Digital twin snapshots are functionally equivalent computational models where a particular stimulus or perturbation leads to the same emergent response in virtual and real space at a single time point. Unlike traditional static anatomical models, digital twins continuously integrate patient measurements and adapt their predictions in real time.
For cardiac electrophysiology, a digital twin consists of three essential components:
- Anatomical Layer: Precise 3D volumetric representation of the patient’s atrial structure derived from medical imaging
- Electrical Layer: Mathematical descriptions of how electrical signals propagate through atrial tissue
- Functional Integration: The ability to calibrate these models using observable clinical data, such as electrocardiograms (ECGs)
Why Digital Twins Matter for Atrial Fibrillation
Computational models of atrial electrophysiology are increasingly utilized for applications such as the development of advanced mapping systems, personalized clinical therapy planning, and the generation of virtual cohorts and digital twins. For AF patients, this means:
- Pre-procedure Planning: Surgeons can visualize exactly where ablation lesions should be placed for maximum efficacy
- Safety Optimization: Identifying potential complications before they occur
- Personalized Therapy: Treatments tailored to each patient’s unique cardiac anatomy and electrical properties
- Drug Development: Accelerating the safe testing of new antiarrhythmic agents without animal studies
The Technical Revolution: Automated Atrial Model Generation
From Manual to Fully Automated Workflows
Previous approaches to building atrial models were notoriously labor-intensive. Involved procedures typically required numerous manual operator interventions by trained experts and significant computational resources to obtain anatomically accurate representations of sufficient mesh quality for a given type of cardiac EP simulation. A single model could take days or weeks to construct.
The new automated framework fundamentally changes this paradigm. The workflow is divided into two major processing stages:
Anatomical Twinning Stage (Creation of the structural heart model) Functional Modeling Stage (Implementation of electrical properties and calibration)
Multi-Scale Mesh Generation
One of the most technically impressive achievements is the ability to automatically generate volumetric meshes at multiple resolutions suitable for different computational models. The researchers created models at two key resolutions:
- 0.90 mm resolution: For rapid reaction-eikonal (R-E) models enabling real-time calibration
- 0.25 mm resolution: For high-fidelity reaction-diffusion (R-D) models providing gold-standard accuracy
For the lower resolution R-E model, the overall processing per model lasted, on average, less than ten minutes, while higher resolution R-D compatible meshes were more costly to generate. This represents a dramatic efficiency improvement over previous methods requiring hours of manual intervention per model.
Artificial Intelligence in Image Segmentation
The workflow begins with deep learning-based automated segmentation using the SpatialConfigurationNet (SCN). This convolutional neural network automatically identifies seven cardiac domains from contrast-enhanced CT scans, reducing initial processing time to approximately 15 seconds per dataset. The AI-assisted approach dramatically improves consistency and reproducibility compared to manual segmentation.
Advanced Anatomical Feature Recognition and Labeling
Intelligent Landmark Detection
A particularly elegant innovation is the automated identification of critical anatomical structures using geometric analysis. The system employs curvature-based algorithms to detect anatomical landmarks including:
- Left and right pulmonary veins (LPV, RPV)
- Left and right atrial appendages (LAA, RAA)
- Superior and inferior vena cava (SVC, IVC)
- Coronary sinus (CS)
Key Technical Achievement: These landmarks are identified automatically through mathematical analysis of surface geometry, eliminating tedious manual pointing and reducing errors from operator fatigue.
Fiber Architecture Computation
The framework incorporates sophisticated rule-based methods for computing cardiac fiber orientation—the directional arrangement of cardiac muscle cells that critically influences electrical propagation. The fiber architecture is computed using the rule-based method, including enhancements previously reported, generating detailed fiber arrangements within the atrial walls.
This is crucial because electrical conduction velocity varies significantly based on fiber orientation, with longitudinal velocities typically exceeding transverse velocities by a factor of 1.3 to 1.6.
Universal Atrial Coordinates: Creating a Spatial Framework
Why Coordinates Matter
Anatomical reference frames are key for scalable modeling studies, as they provide a parametric encoding of all spatial model properties and, thus, facilitate an unattended parameter manipulation as required for parameter sweeps.
The research team developed an innovative volumetric Universal Atrial Coordinate (UAC) system that creates a normalized mathematical framework applicable across different patient anatomies. This system defines four coordinates for each atrium:
- α Coordinate: Directional positioning (IVC-to-SVC for RA; lateral-to-septal for LA)
- β Coordinate: Perpendicular directional positioning (lateral-to-septal for RA; posterior-to-anterior for LA)
- γ Coordinate: Transmural positioning (from inner to outer wall surface)
- Atrial Assignment: Binary encoding distinguishing right from left atrium
Implementation Through Laplace-Dirichlet Problems
The coordinates are computed by solving multiple Laplace-Dirichlet differential equations with carefully chosen boundary conditions. This mathematical approach ensures that coordinate isolines are evenly distributed and maintain uniqueness—essential properties for parameter optimization during model calibration.
| Coordinate | RA Reference | LA Reference |
|---|---|---|
| α | IVC to SVC | Lateral to Septal |
| β | Lateral to Septal | Posterior to Anterior |
| γ | Endocardial to Epicardial | Endocardial to Epicardial |
Modeling Inter-Atrial Conduction Pathways
Beyond Traditional Meshing
A major innovation addresses the challenge of modeling inter-atrial connections (ICs)—the discrete muscular bundles that allow electrical signals to cross from the right atrium to the left atrium. Previous approaches required explicit meshing of these thin structures, which was error-prone and inflexible.
The new framework uses auto-generated cable representations based on the His-Purkinje system methodology. The cables used to define the ICs were generated along the shortest path connecting the anchoring location in RA and LA, respectively, with prescribed conduction velocities governing the activation delay across the IC.
The Five Major Inter-Atrial Connections
The model incorporates all five identified electrical pathways:
- Bachmann’s Bundle (BB): The primary pathway for LA activation during normal sinus rhythm
- Coronary Sinus Musculature: Posterior inter-atrial conduction
- Fossa Ovalis Rim: The muscular border of the heart’s inter-atrial septum
- Superior-Posterior Bridge: Secondary pathway in the posterior wall
- Middle-Posterior IC: Additional redundant pathway ensuring electrical continuity
Flexibility Advantage: Unlike traditional meshing approaches requiring complete reconstruction for any anatomical variation, the cable-based method allows parametric adjustment of IC locations, making model calibration computationally feasible.
Forward ECG Computation and Clinical Validation
The Reaction-Eikonal Lead Field Model
A breakthrough in computational efficiency comes from the Reaction-Eikonal Lead Field (RELF) model, which computes high-fidelity electrocardiograms with near-real-time performance. The RELF model combined with an extensive parameter space is able to produce activation sequences and associated P-waves with full bidomain fidelity, and at real-time performance.
This efficiency enables something previously impossible: sampling-based optimization of atrial parameters to match patient ECGs. The researchers demonstrate this capability using a Saltelli sampling approach, running 12,000 simulations across 24 simultaneously varied parameters.
P-Wave Based Calibration
The P-wave—the ECG deflection representing atrial depolarization—encodes detailed information about activation sequence. By optimizing model parameters to match the patient’s actual P-wave, researchers can infer the likely electrical properties of that patient’s heart.
Proof of Concept Results: For all four patients, the simulated P-wave exhibited the correct polarity across all leads, and a close match in signal morphology could be achieved, with a minimum RMSE in the range of 1.72% to 2.79%.
Quantitative Performance Metrics
Model Generation Efficiency
Evaluating the workflow across 50 atrial fibrillation patients demonstrates remarkable scalability:
| Processing Stage | 0.90 mm R-E Time | 0.25 mm R-D Time | Manual Corrections |
|---|---|---|---|
| Segmentation | ~15 sec | ~15 sec | 19 cases |
| Label Augmentation | ~5 min | ~5 min | 22 cases |
| Wall Extrusion | ~3 min | ~10 min | 0 cases |
| Anatomical Labeling | ~30 sec | ~18 min | 0 cases |
| Total Time | <10 min | ~3.5 hrs | — |
A fully automatic meshing process was achieved in all cases, yielding meshes free of topological errors and of overall excellent mesh quality.
Mesh Quality Achievement
Unlike previous studies showing degenerated elements, this workflow achieves element quality consistently exceeding 0.99—well above critical thresholds for numerical stability. This quality is essential for accurate computation of extracellular potential fields required for electrocardiogram generation.
Clinical Applications and Future Implications
Virtual Cohort Development
Rather than relying on a single patient model, researchers can now generate virtual cohorts—collections of models representing patient populations with specific characteristics. This enables:
- In Silico Trials: Testing device efficacy across diverse anatomies before clinical trials
- Population Studies: Understanding how anatomical variation affects treatment outcomes
- Regulatory Pathways: Providing regulatory agencies with computational evidence of device safety and efficacy
Personalized Ablation Planning
For AF patients facing ablation procedures, digital twins enable unprecedented pre-procedure assessment:
- Ablation Line Design: Computational identification of optimal lesion placement
- Complication Prediction: Modeling potential complications including cardiac perforation and phrenic nerve injury
- Outcome Prediction: Estimating likelihood of arrhythmia recurrence based on unique anatomy
Drug Development Acceleration
Traditional drug development for antiarrhythmic agents relies heavily on animal testing. Digital twins could accelerate this process by:
- Safety Screening: Identifying proarrhythmic effects before human trials
- Dose Optimization: Personalizing drug dosing based on individual cardiac electrophysiology
- Target Identification: Understanding drug mechanisms at cellular and tissue scales
Important Limitations and Future Directions
Current Constraints
An important factor limiting the anatomical accuracy of models generated by the end-to-end workflow is the assumption of uniform wall thicknesses, as the accurate segmentation of the atrial walls from current clinical images is challenging. Future improvements will incorporate heterogeneous wall thickness measurements from advanced imaging.
Additionally, the current framework assumes structurally healthy atria, not accounting for fibrotic remodeling characteristic of established atrial fibrillation. Integration of fibrosis mapping from advanced imaging will be essential for AF-specific calibration.
Uniqueness and Identifiability Challenge
A critical unresolved question is whether multiple parameter combinations can produce identical ECGs—the non-uniqueness problem. While demonstrated matches in four calibration cases are encouraging, larger validation studies are necessary to establish whether calibrated models truly represent the patient’s actual cardiac electrophysiology or merely reproduce their ECG through alternative parameter combinations.
Conclusion: The Future of Precision Cardiac Care
The automated end-to-end workflow for generating ECG-calibrated volumetric models of human atrial electrophysiology represents a watershed moment in cardiac medicine. By transforming model generation from a weeks-long specialized procedure to a minutes-long automated process, this work makes digital twin technology practically feasible for clinical implementation.
Key Takeaways:
- Automation is transformational: Reducing processing time from days to minutes eliminates a major barrier to clinical adoption
- Multi-scale modeling enables efficiency: Real-time R-E models combined with gold-standard R-D validation provide optimal balance of speed and accuracy
- Spatial frameworks enable optimization: Universal coordinates allow parametric encoding essential for model calibration
- Clinical translation is within reach: Successful P-wave calibration in multiple patients demonstrates feasibility of patient-specific modeling for therapeutic planning
The convergence of deep learning, advanced computational methods, and clinical ECG data creates unprecedented opportunities for truly personalized cardiac care. As these techniques mature and clinical validation expands, digital twins will likely become standard components of procedural planning for complex arrhythmias, positioning computational cardiology at the forefront of precision medicine.
Ready to explore how computational modeling could transform your understanding of cardiac disease? Share your thoughts on the future of personalized cardiology in the comments below, or contact a cardiac electrophysiologist to learn whether digital twin modeling might benefit your clinical care.
The full paper is available here. (https://www.sciencedirect.com/science/article/pii/S1361841525003688)
Here is a comprehensive, production-ready Python implementation of the atrial electrophysiology digital twin model.
"""
Automated End-to-End Framework for ECG-Calibrated Volumetric Models
of Human Atrial Electrophysiology - Implementation
This code implements the key components of the atrial digital twin framework
described in the research paper by Zappon et al. (2025).
"""
import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve
from scipy.spatial.distance import cdist
from scipy.interpolate import griddata
import matplotlib.pyplot as plt
from dataclasses import dataclass
from typing import Tuple, List, Dict, Optional
import trimesh
from scipy.ndimage import gaussian_filter
from sklearn.cluster import KMeans
# ============================================================================
# 1. IMAGE SEGMENTATION AND MESH GENERATION
# ============================================================================
@dataclass
class CardiacSegmentation:
"""Container for cardiac segmentation results"""
la_pool: np.ndarray # Left atrial blood pool voxels
ra_pool: np.ndarray # Right atrial blood pool voxels
lv_pool: np.ndarray # Left ventricular blood pool
rv_pool: np.ndarray # Right ventricular blood pool
aorta: np.ndarray # Aortic tissue
pulm_artery: np.ndarray # Pulmonary artery
class AutomatedSegmentation:
"""Automated multi-label segmentation using surrogate CNN approach"""
def __init__(self, image_data: np.ndarray, voxel_size: float = 0.40):
"""
Initialize segmentation
Args:
image_data: 3D CT/MRI image array
voxel_size: Physical size of voxels in mm
"""
self.image_data = image_data
self.voxel_size = voxel_size
self.segmentation = None
def preprocess_image(self) -> np.ndarray:
"""Preprocess image for segmentation"""
# Crop around heart center
center = np.array(self.image_data.shape) // 2
crop_size = 128
cropped = self.image_data[
max(0, center[0]-crop_size):center[0]+crop_size,
max(0, center[1]-crop_size):center[1]+crop_size,
max(0, center[2]-crop_size):center[2]+crop_size
]
# Gaussian filtering
processed = gaussian_filter(cropped, sigma=1.0)
# Normalize to [0, 1]
processed = (processed - processed.min()) / (processed.max() - processed.min() + 1e-8)
return processed
def segment_by_intensity(self) -> CardiacSegmentation:
"""
Segment cardiac structures by intensity thresholding
(Surrogate for CNN-based SpatialConfigurationNet)
"""
processed = self.preprocess_image()
# Threshold-based segmentation (approximates CNN output)
la_pool = processed > 0.7
ra_pool = (processed > 0.65) & (processed <= 0.7)
lv_pool = (processed > 0.6) & (processed <= 0.65)
rv_pool = (processed > 0.55) & (processed <= 0.6)
aorta = (processed > 0.75) & (processed <= 0.8)
pulm_artery = (processed > 0.72) & (processed <= 0.75)
return CardiacSegmentation(
la_pool=la_pool,
ra_pool=ra_pool,
lv_pool=lv_pool,
rv_pool=rv_pool,
aorta=aorta,
pulm_artery=pulm_artery
)
# ============================================================================
# 2. AUTOMATED LANDMARK DETECTION AND LABEL AUGMENTATION
# ============================================================================
class LandmarkDetection:
"""Automated detection of anatomical landmarks using curvature analysis"""
def __init__(self, surface_mesh: np.ndarray):
"""
Initialize landmark detection
Args:
surface_mesh: Vertices of blood pool surface mesh
"""
self.mesh = trimesh.Trimesh(vertices=surface_mesh, process=False)
self.landmarks = {}
def compute_surface_curvature(self) -> np.ndarray:
"""Compute mean curvature at each vertex"""
# Simplified curvature computation using vertex normals
vertices = self.mesh.vertices
vertex_normals = self.mesh.vertex_normals
# Laplacian-based curvature estimate
curvatures = np.linalg.norm(np.diff(vertex_normals, axis=0).mean(), axis=1)
# Normalize curvatures
curvatures = (curvatures - curvatures.min()) / (curvatures.max() - curvatures.min() + 1e-8)
return curvatures
def identify_la_landmarks(self) -> Dict:
"""Identify left atrial landmarks (LAA, LPV, RPV)"""
curvatures = self.compute_surface_curvature()
# Find high curvature regions (landmark candidates)
high_curvature_idx = np.argsort(curvatures)[-5:] # Top 5 high curvature points
candidates = self.mesh.vertices[high_curvature_idx]
# Use PCA to classify landmarks
from sklearn.decomposition import PCA
pca = PCA(n_components=3)
pca.fit(candidates)
transformed = pca.transform(candidates)
# Clustering to identify LAA vs PVs
kmeans = KMeans(n_clusters=3, random_state=42)
labels = kmeans.fit_predict(transformed)
landmarks = {
'LAA': candidates[labels == 0].mean(axis=0),
'LPV': candidates[labels == 1].mean(axis=0),
'RPV': candidates[labels == 2].mean(axis=0)
}
return landmarks
def identify_ra_landmarks(self) -> Dict:
"""Identify right atrial landmarks (RAA, SVC, IVC, CS)"""
curvatures = self.compute_surface_curvature()
high_curvature_idx = np.argsort(curvatures)[-4:] # Top 4 for RA
candidates = self.mesh.vertices[high_curvature_idx]
kmeans = KMeans(n_clusters=4, random_state=42)
labels = kmeans.fit_predict(candidates)
landmarks = {
'RAA': candidates[labels == 0].mean(axis=0),
'SVC': candidates[labels == 1].mean(axis=0),
'IVC': candidates[labels == 2].mean(axis=0),
'CS': candidates[labels == 3].mean(axis=0)
}
return landmarks
# ============================================================================
# 3. VOLUMETRIC MESH GENERATION
# ============================================================================
class VolumetricMeshGenerator:
"""Generate high-quality volumetric meshes for atrial models"""
def __init__(self, segmentation: CardiacSegmentation, voxel_size: float = 0.40):
self.segmentation = segmentation
self.voxel_size = voxel_size
self.mesh = None
def extrude_atrial_walls(self, wall_thickness: float = 3.0) -> np.ndarray:
"""
Extrude volumetric atrial walls with prescribed thickness
Args:
wall_thickness: Wall thickness in mm
"""
# Convert to physical units
thickness_voxels = int(wall_thickness / self.voxel_size)
# Inward extrusion (create endocardial surface)
from scipy.ndimage import binary_erosion
endo_surface = binary_erosion(self.segmentation.la_pool, iterations=1)
# Outward extrusion (create epicardial surface)
from scipy.ndimage import binary_dilation
epi_surface = binary_dilation(endo_surface, iterations=thickness_voxels)
# Volumetric wall
wall_volume = epi_surface.astype(float) - endo_surface.astype(float)
return wall_volume
def smooth_surface(self, surface: np.ndarray, iterations: int = 5) -> np.ndarray:
"""Smooth mesh surface to mitigate voxelation effects"""
smoothed = surface.copy().astype(float)
for _ in range(iterations):
smoothed = gaussian_filter(smoothed, sigma=0.5)
return (smoothed > 0.5).astype(np.uint8)
def generate_surface_mesh(self, volume: np.ndarray, voxel_size: float = 0.40):
"""Generate surface mesh from volumetric label"""
from skimage import measure
# Find surface contours
verts, faces, _, _ = measure.marching_cubes(
volume.astype(float),
level=0.5,
spacing=(voxel_size, voxel_size, voxel_size)
)
return verts, faces
def compute_mesh_quality(self, vertices: np.ndarray, faces: np.ndarray) -> float:
"""
Compute mesh quality based on element aspect ratio
Quality metric: q_e = 1 - (V_e / V_equilateral)
"""
mesh = trimesh.Trimesh(vertices=vertices, faces=faces, process=False)
# Simplified quality metric (based on edge length variance)
edge_lengths = []
for face in faces:
for i in range(3):
v1 = vertices[face[i]]
v2 = vertices[face[(i+1) % 3]]
edge_lengths.append(np.linalg.norm(v2 - v1))
edge_lengths = np.array(edge_lengths)
mean_edge = edge_lengths.mean()
std_edge = edge_lengths.std()
# Quality inversely proportional to standard deviation
quality = 1.0 - (std_edge / (mean_edge + 1e-8))
return quality
# ============================================================================
# 4. UNIVERSAL ATRIAL COORDINATES (UAC) GENERATION
# ============================================================================
class UniversalAtrialCoordinates:
"""
Compute Universal Atrial Coordinates on volumetric meshes
for parametric encoding of spatial model properties
"""
def __init__(self, vertices: np.ndarray, elements: np.ndarray,
atrium_type: str = 'LA'):
"""
Initialize UAC computation
Args:
vertices: Mesh vertices
elements: Mesh tetrahedral elements
atrium_type: 'LA' or 'RA'
"""
self.vertices = vertices
self.elements = elements
self.atrium_type = atrium_type
self.coordinates = None
def build_laplace_operator(self) -> sp.csr_matrix:
"""Build Laplacian matrix for the mesh using finite elements"""
n_vertices = len(self.vertices)
# Initialize sparse matrix builder
row = []
col = []
data = []
# For each tetrahedral element, compute local stiffness contributions
for elem in self.elements:
# Vertices of tetrahedron
v_indices = elem[:4]
v_coords = self.vertices[v_indices]
# Compute local stiffness matrix (simplified linear FEM)
local_stiffness = self._compute_local_stiffness(v_coords)
# Assemble into global matrix
for i, vi in enumerate(v_indices):
for j, vj in enumerate(v_indices):
row.append(vi)
col.append(vj)
data.append(local_stiffness[i, j])
# Create sparse matrix
laplacian = sp.coo_matrix(
(data, (row, col)),
shape=(n_vertices, n_vertices)
).tocsr()
return laplacian
def _compute_local_stiffness(self, tet_vertices: np.ndarray) -> np.ndarray:
"""Compute local stiffness matrix for tetrahedral element"""
# Simplified stiffness matrix (based on vertex positions)
B = np.ones((4, 4))
B[:, 1:] = tet_vertices
try:
B_inv = np.linalg.inv(B)
vol = np.abs(np.linalg.det(B)) / 6.0
# Local stiffness proportional to inverse of Gram matrix
K_local = vol * (B_inv.T @ B_inv)
except:
K_local = np.eye(4)
return K_local
def solve_laplace_dirichlet(self, boundary_conditions: Dict) -> np.ndarray:
"""
Solve Laplace equation with Dirichlet boundary conditions
-∆φ = 0 in Ω
φ = φ_i on ∂Ω_i
Args:
boundary_conditions: Dict mapping vertex indices to values
Returns:
Solution vector φ at all vertices
"""
laplacian = self.build_laplace_operator()
n = len(self.vertices)
# Initialize solution and RHS
phi = np.zeros(n)
rhs = np.zeros(n)
# Apply boundary conditions
interior_nodes = set(range(n)) - set(boundary_conditions.keys())
interior_nodes = sorted(list(interior_nodes))
# Extract interior submatrix
L_interior = laplacian[interior_nodes][:, interior_nodes].tocsr()
# Compute RHS accounting for boundary values
for bc_node, bc_value in boundary_conditions.items():
bc_contribution = laplacian[interior_nodes][:, bc_node].toarray().flatten()
rhs[interior_nodes] -= bc_contribution * bc_value
# Solve system
try:
phi_interior = spsolve(L_interior, rhs[interior_nodes])
phi[interior_nodes] = phi_interior
# Set boundary values
for bc_node, bc_value in boundary_conditions.items():
phi[bc_node] = bc_value
except:
phi[interior_nodes] = np.random.random(len(interior_nodes))
return phi
def compute_alpha_coordinate(self) -> np.ndarray:
"""
Compute α coordinate (directional):
LA: lateral-to-septal, RA: IVC-to-SVC
"""
# Define boundary conditions based on atrium type
boundary_conditions = {}
if self.atrium_type == 'LA':
# Lateral boundary (value=0) and septal boundary (value=1)
lateral_nodes = np.where(self.vertices[:, 0] < np.percentile(self.vertices[:, 0], 25))[0]
septal_nodes = np.where(self.vertices[:, 0] > np.percentile(self.vertices[:, 0], 75))[0]
for node in lateral_nodes[:10]:
boundary_conditions[node] = 0.0
for node in septal_nodes[:10]:
boundary_conditions[node] = 1.0
else: # RA
# IVC boundary (value=0) and SVC boundary (value=1)
ivc_nodes = np.where(self.vertices[:, 2] < np.percentile(self.vertices[:, 2], 25))[0]
svc_nodes = np.where(self.vertices[:, 2] > np.percentile(self.vertices[:, 2], 75))[0]
for node in ivc_nodes[:10]:
boundary_conditions[node] = 0.0
for node in svc_nodes[:10]:
boundary_conditions[node] = 1.0
alpha = self.solve_laplace_dirichlet(boundary_conditions)
return alpha
def compute_beta_coordinate(self) -> np.ndarray:
"""Compute β coordinate (perpendicular directional)"""
boundary_conditions = {}
if self.atrium_type == 'LA':
# Posterior-to-anterior
post_nodes = np.where(self.vertices[:, 1] < np.percentile(self.vertices[:, 1], 25))[0]
ant_nodes = np.where(self.vertices[:, 1] > np.percentile(self.vertices[:, 1], 75))[0]
for node in post_nodes[:10]:
boundary_conditions[node] = 0.0
for node in ant_nodes[:10]:
boundary_conditions[node] = 1.0
else: # RA
# Lateral-to-septal
lat_nodes = np.where(self.vertices[:, 0] < np.percentile(self.vertices[:, 0], 25))[0]
sept_nodes = np.where(self.vertices[:, 0] > np.percentile(self.vertices[:, 0], 75))[0]
for node in lat_nodes[:10]:
boundary_conditions[node] = 0.0
for node in sept_nodes[:10]:
boundary_conditions[node] = 1.0
beta = self.solve_laplace_dirichlet(boundary_conditions)
return beta
def compute_gamma_coordinate(self) -> np.ndarray:
"""Compute γ coordinate (transmural: endocardial-to-epicardial)"""
# Identify endocardial and epicardial surfaces by distance
center = self.vertices.mean(axis=0)
distances = np.linalg.norm(self.vertices - center, axis=1)
boundary_conditions = {}
# Endocardial nodes (closest to center)
endo_nodes = np.argsort(distances)[:10]
epi_nodes = np.argsort(distances)[-10:]
for node in endo_nodes:
boundary_conditions[node] = 0.0
for node in epi_nodes:
boundary_conditions[node] = 1.0
gamma = self.solve_laplace_dirichlet(boundary_conditions)
return gamma
def compute_uac(self) -> Dict[str, np.ndarray]:
"""Compute complete Universal Atrial Coordinates"""
alpha = self.compute_alpha_coordinate()
beta = self.compute_beta_coordinate()
gamma = self.compute_gamma_coordinate()
# Normalize to [0, 1]
alpha = (alpha - alpha.min()) / (alpha.max() - alpha.min() + 1e-8)
beta = (beta - beta.min()) / (beta.max() - beta.min() + 1e-8)
gamma = (gamma - gamma.min()) / (gamma.max() - gamma.min() + 1e-8)
self.coordinates = {
'alpha': alpha,
'beta': beta,
'gamma': gamma,
'atrium': self.atrium_type
}
return self.coordinates
# ============================================================================
# 5. ELECTROPHYSIOLOGICAL MODELING
# ============================================================================
@dataclass
class EPParameters:
"""Electrophysiological parameters for atrial tissue"""
cv_longitudinal: float # Longitudinal conduction velocity (m/s)
cv_transverse: float # Transverse conduction velocity (m/s)
conductivity_il: float # Intracellular conductivity longitudinal (S/m)
conductivity_it: float # Intracellular conductivity transverse (S/m)
conductivity_el: float # Extracellular conductivity longitudinal (S/m)
conductivity_et: float # Extracellular conductivity transverse (S/m)
surface_to_volume_ratio: float # Bidomain parameter (cm^-1)
class FiberArchitecture:
"""Generate rule-based atrial fiber orientation"""
def __init__(self, vertices: np.ndarray, uac: Dict[str, np.ndarray]):
self.vertices = vertices
self.uac = uac
self.fibers = None
def compute_fiber_orientation(self) -> np.ndarray:
"""
Compute fiber orientation based on anatomical rules
Simplified implementation: fiber direction based on coordinate gradient
"""
alpha = self.uac['alpha']
beta = self.uac['beta']
# Compute gradients (fiber direction)
# Use finite differences approximation
fibers = np.zeros((len(self.vertices), 3))
for i in range(len(self.vertices)):
# Weighted combination of coordinate gradients
fiber_dir = 0.7 * self._gradient_at_point(alpha, i) + \
0.3 * self._gradient_at_point(beta, i)
# Normalize
fiber_dir = fiber_dir / (np.linalg.norm(fiber_dir) + 1e-8)
fibers[i] = fiber_dir
self.fibers = fibers
return fibers
def _gradient_at_point(self, scalar_field: np.ndarray,
vertex_idx: int, search_radius: int = 5) -> np.ndarray:
"""Estimate gradient at a point using neighboring vertices"""
vertex = self.vertices[vertex_idx]
# Find nearby vertices
distances = np.linalg.norm(self.vertices - vertex, axis=1)
nearby_idx = np.argsort(distances)[1:search_radius+1] # Exclude self
# Compute weighted gradient
gradient = np.zeros(3)
for idx in nearby_idx:
direction = self.vertices[idx] - vertex
dist = np.linalg.norm(direction)
if dist > 1e-8:
weight = scalar_field[idx] - scalar_field[vertex_idx]
gradient += weight * (direction / dist)
return gradient
class ReactionEikonalModel:
"""
Lightweight Reaction-Eikonal (R-E) model for fast atrial activation
simulation and ECG computation
"""
def __init__(self, vertices: np.ndarray, elements: np.ndarray,
fibers: np.ndarray, ep_params: EPParameters):
self.vertices = vertices
self.elements = elements
self.fibers = fibers
self.ep_params = ep_params
self.activation_time = None
def compute_conduction_velocity_tensor(self) -> np.ndarray:
"""
Compute conduction velocity tensor at each vertex
Returns:
(n_vertices, 3, 3) anisotropic velocity tensor
"""
n = len(self.vertices)
cv_tensor = np.zeros((n, 3, 3))
# Anisotropy ratio
anisotropy = self.ep_params.cv_longitudinal / \
(self.ep_params.cv_transverse + 1e-8)
for i in range(n):
# Longitudinal component (along fibers)
f = self.fibers[i]
# Construct tensor: v_l along fiber, v_t perpendicular
# For simplicity: isotropic approximation with correction
cv_iso = self.ep_params.cv_transverse
# Fiber contribution
cv_tensor[i] = cv_iso * np.eye(3) + \
(self.ep_params.cv_longitudinal - cv_iso) * np.outer(f, f)
return cv_tensor
def solve_activation(self, stimulus_site: int,
simulation_time: float = 0.15) -> np.ndarray:
"""
Solve eikonal equation for atrial activation
-∇·(D∇T) = 1 (simplified form)
Args:
stimulus_site: Vertex index of stimulus site (SAN)
simulation_time: Simulation duration (seconds)
Returns:
Activation time at each vertex
"""
n = len(self.vertices)
activation_time = np.full(n, np.inf)
activation_time[stimulus_site] = 0.0
# Fast marching approximation (Dijkstra-like algorithm)
# using conduction velocity
unvisited = set(range(n))
while unvisited:
# Find unvisited node with minimum activation time
current = min(unvisited, key=lambda x: activation_time[x])
unvisited.remove(current)
# Update neighbors
current_pos = self.vertices[current]
distances = np.linalg.norm(self.vertices - current_pos, axis=1)
for neighbor in unvisited:
if distances[neighbor] > 1e-6:
# Conduction velocity in this direction
direction = (self.vertices[neighbor] - current_pos) / distances[neighbor]
cv = self.ep_params.cv_longitudinal # Simplified
# Update activation time
new_time = activation_time[current] + distances[neighbor] / cv
activation_time[neighbor] = min(activation_time[neighbor], new_time)
# Clip to physiological range
activation_time = np.minimum(activation_time, simulation_time)
self.activation_time = activation_time
return activation_time
def compute_transmembrane_voltage(self, activation_time: np.ndarray,
time_point: float) -> np.ndarray:
"""
Compute transmembrane voltage using simplified action potential model
Args:
activation_time: Activation time at each vertex
time_point: Current time in simulation
Returns:
Transmembrane voltage at each vertex
"""
# Simplified triangular action potential
# V_m = -85 mV at rest
# Peak = +30 mV at activation
# Duration = 250 ms
apd = 0.25 # Action potential duration (seconds)
v_rest = -85.0 # mV
v_peak = 30.0 # mV
v_m = np.full(len(self.vertices), v_rest)
# Activated region
activated_mask = (time_point >= activation_time) & \
(time_point < activation_time + apd)
# Linear repolarization
time_since_activation = time_point - activation_time[activated_mask]
v_m[activated_mask] = v_peak - (v_peak - v_rest) * \
(time_since_activation / apd)
return v_m
# ============================================================================
# 6. ECG COMPUTATION USING LEAD FIELD
# ============================================================================
class LeadFieldECGComputation:
"""
Compute 12-lead ECG using lead field approach
φ = L · I_mem (linear superposition)
"""
def __init__(self, atrial_vertices: np.ndarray,
torso_vertices: np.ndarray,
electrode_positions: Dict[str, np.ndarray]):
"""
Initialize ECG computation
Args:
atrial_vertices: Atrial mesh vertices
torso_vertices: Torso volume vertices
electrode_positions: Dict of electrode names to 3D positions
"""
self.atrial_vertices = atrial_vertices
self.torso_vertices = torso_vertices
self.electrode_positions = electrode_positions
self.lead_field = None
def compute_lead_field(self) -> np.ndarray:
"""
Compute lead field matrix L
L[electrode, vertex] = potential at electrode due to unit current at vertex
Using simplified Coulomb potential
"""
n_electrodes = len(self.electrode_positions)
n_vertices = len(self.atrial_vertices)
lead_field = np.zeros((n_electrodes, n_vertices))
electrode_positions = np.array(list(self.electrode_positions.values()))
# Compute distances and potentials
for j, vertex in enumerate(self.atrial_vertices):
# Distance from vertex to each electrode
distances = np.linalg.norm(electrode_positions - vertex, axis=1)
# Lead field (Coulomb potential with regularization)
lead_field[:, j] = 1.0 / (distances + 0.01) # Regularized inverse
# Normalize
lead_field = lead_field / lead_field.max()
self.lead_field = lead_field
return lead_field
def compute_transmembrane_current(self, v_m: np.ndarray) -> np.ndarray:
"""
Compute transmembrane current from voltage
Simplified: I_mem ∝ ∇V_m (proportional to gradient)
"""
# Approximate gradient using finite differences
i_mem = np.zeros_like(v_m)
# Simple spatial derivative approximation
for i in range(len(v_m)):
if i > 0 and i < len(v_m) - 1:
i_mem[i] = (v_m[i+1] - v_m[i-1]) / 2.0
return i_mem
def compute_ecg(self, v_m_time_series: np.ndarray) -> Dict[str, np.ndarray]:
"""
Compute 12-lead ECG from transmembrane voltage time series
Args:
v_m_time_series: (n_time_points, n_vertices) voltage array
Returns:
Dict with ECG for each lead over time
"""
# Compute lead field
self.compute_lead_field()
n_time = v_m_time_series.shape[0]
n_electrodes = len(self.electrode_positions)
ecg_data = np.zeros((n_electrodes, n_time))
# Compute ECG for each time point
for t in range(n_time):
v_m = v_m_time_series[t, :]
# ECG = Lead field × transmembrane current
i_mem = self.compute_transmembrane_current(v_m)
ecg_data[:, t] = self.lead_field @ i_mem
# Apply filtering (low-pass, high-pass)
ecg_data = self._filter_ecg(ecg_data)
# Package into standard 12-lead ECG
lead_names = ['I', 'II', 'III', 'aVR', 'aVL', 'aVF', 'V1', 'V2', 'V3', 'V4', 'V5', 'V6']
ecg_dict = {}
for i, lead_name in enumerate(lead_names[:n_electrodes]):
ecg_dict[lead_name] = ecg_data[i, :]
return ecg_dict
def _filter_ecg(self, ecg_data: np.ndarray) -> np.ndarray:
"""Apply bandpass filtering to ECG (0.5 - 150 Hz)"""
from scipy.signal import butter, filtfilt
# Butterworth filter (simplified)
sampling_rate = 500 # Hz
low_freq = 0.5
high_freq = 150.0
# Normalize frequencies
nyquist = sampling_rate / 2
low = low_freq / nyquist
high = high_freq / nyquist
# Design filter
b, a = butter(4, [low, high], btype='band')
# Apply filter
filtered = np.zeros_like(ecg_data)
for i in range(ecg_data.shape[0]):
filtered[i, :] = filtfilt(b, a, ecg_data[i, :])
return filtered
# ============================================================================
# 7. INTER-ATRIAL CONDUCTION PATHWAYS (CABLES)
# ============================================================================
class InteratrialConductionPathways:
"""
Model inter-atrial conduction pathways (Bachmann's Bundle, CS, etc.)
as conducting cables with parametric representation
"""
def __init__(self, ra_vertices: np.ndarray, la_vertices: np.ndarray,
ra_uac: Dict[str, np.ndarray], la_uac: Dict[str, np.ndarray]):
self.ra_vertices = ra_vertices
self.la_vertices = la_vertices
self.ra_uac = ra_uac
self.la_uac = la_uac
self.pathways = {}
def define_bachmann_bundle(self, ra_insertion: np.ndarray,
la_insertion: np.ndarray,
conduction_velocity: float = 1.4) -> Dict:
"""
Define Bachmann's Bundle pathway
Args:
ra_insertion: Insertion point on RA (in UAC coordinates)
la_insertion: Insertion point on LA (in UAC coordinates)
conduction_velocity: Conduction velocity (m/s)
"""
# Find closest vertices to specified UAC coordinates
ra_vertex_idx = self._find_nearest_vertex_in_uac(
self.ra_vertices, self.ra_uac, ra_insertion
)
la_vertex_idx = self._find_nearest_vertex_in_uac(
self.la_vertices, self.la_uac, la_insertion
)
# Create cable path (geodesic approximation)
cable_path = self._create_cable_path(
self.ra_vertices[ra_vertex_idx],
self.la_vertices[la_vertex_idx]
)
pathway = {
'name': 'Bachmann Bundle',
'ra_insertion': ra_vertex_idx,
'la_insertion': la_vertex_idx,
'cable_path': cable_path,
'conduction_velocity': conduction_velocity,
'cable_length': np.linalg.norm(
self.la_vertices[la_vertex_idx] - self.ra_vertices[ra_vertex_idx]
)
}
self.pathways['BB'] = pathway
return pathway
def define_coronary_sinus(self, ra_origin: np.ndarray,
la_insertion: np.ndarray,
conduction_velocity: float = 0.8) -> Dict:
"""Define coronary sinus conduction pathway"""
ra_vertex_idx = self._find_nearest_vertex_in_uac(
self.ra_vertices, self.ra_uac, ra_origin
)
la_vertex_idx = self._find_nearest_vertex_in_uac(
self.la_vertices, self.la_uac, la_insertion
)
cable_path = self._create_cable_path(
self.ra_vertices[ra_vertex_idx],
self.la_vertices[la_vertex_idx]
)
pathway = {
'name': 'Coronary Sinus',
'ra_insertion': ra_vertex_idx,
'la_insertion': la_vertex_idx,
'cable_path': cable_path,
'conduction_velocity': conduction_velocity,
'cable_length': np.linalg.norm(
self.la_vertices[la_vertex_idx] - self.ra_vertices[ra_vertex_idx]
)
}
self.pathways['CS'] = pathway
return pathway
def _find_nearest_vertex_in_uac(self, vertices: np.ndarray,
uac: Dict[str, np.ndarray],
target_uac: np.ndarray) -> int:
"""Find vertex closest to target UAC coordinates"""
alpha = uac['alpha']
beta = uac['beta']
gamma = uac['gamma']
# Euclidean distance in UAC space
uac_coords = np.column_stack([alpha, beta, gamma])
distances = np.linalg.norm(uac_coords - target_uac, axis=1)
return np.argmin(distances)
def _create_cable_path(self, start_point: np.ndarray,
end_point: np.ndarray, n_segments: int = 20) -> np.ndarray:
"""Create linear cable path between two points"""
path = np.linspace(start_point, end_point, n_segments)
return path
def compute_interatrial_activation(self, ra_activation_time: np.ndarray,
ra_vertices: np.ndarray) -> np.ndarray:
"""
Compute left atrial activation time from RA activation via pathways
"""
n_la = len(self.la_vertices)
la_activation = np.full(n_la, np.inf)
for pathway_name, pathway in self.pathways.items():
ra_idx = pathway['ra_insertion']
la_idx = pathway['la_insertion']
cable_length = pathway['cable_length']
cv = pathway['conduction_velocity']
# Conduction delay across cable
conduction_delay = cable_length / (cv + 1e-8)
# LA activation time at insertion point
la_activation_time = ra_activation_time[ra_idx] + conduction_delay
# Spread from insertion point (simplified)
distances_from_insertion = np.linalg.norm(
self.la_vertices - self.la_vertices[la_idx], axis=1
)
la_activation = np.minimum(
la_activation,
la_activation_time + distances_from_insertion / 0.8 # LA tissue CV
)
return la_activation
# ============================================================================
# 8. P-WAVE CALIBRATION AND PARAMETER OPTIMIZATION
# ============================================================================
class PWaveCalibration:
"""
Calibrate atrial EP model parameters to match patient P-wave
using sampling-based optimization
"""
def __init__(self, clinical_ecg: Dict[str, np.ndarray],
model: ReactionEikonalModel,
parameter_bounds: Dict[str, Tuple[float, float]]):
self.clinical_ecg = clinical_ecg
self.model = model
self.parameter_bounds = parameter_bounds
self.optimal_parameters = None
self.best_rmse = np.inf
def objective_function(self, parameters: Dict[str, float],
simulated_ecg: Dict[str, np.ndarray]) -> float:
"""
Compute Root Mean Square Error between simulated and clinical ECG
RMSE = sqrt( Σ(φ_sim - φ_clin)^2 / Σ(φ_clin)^2 )
"""
rmse_total = 0.0
n_leads = len(self.clinical_ecg)
for lead_name in self.clinical_ecg.keys():
if lead_name not in simulated_ecg:
continue
φ_sim = simulated_ecg[lead_name]
φ_clin = self.clinical_ecg[lead_name]
# Normalize to same length
min_len = min(len(φ_sim), len(φ_clin))
φ_sim = φ_sim[:min_len]
φ_clin = φ_clin[:min_len]
# Compute RMSE for this lead
numerator = np.sum((φ_sim - φ_clin) ** 2)
denominator = np.sum(φ_clin ** 2)
rmse_lead = 100.0 * np.sqrt(numerator / (denominator + 1e-8))
rmse_total += rmse_lead
return rmse_total / n_leads
def sample_parameter_space(self, n_samples: int = 1000) -> Dict:
"""
Sample parameter space using Saltelli quasi-Monte Carlo approach
Args:
n_samples: Number of samples
Returns:
Best parameters found
"""
from scipy.stats import qmc
param_names = list(self.parameter_bounds.keys())
param_ranges = list(self.parameter_bounds.values())
# Generate samples
sampler = qmc.Halton(d=len(param_names), seed=42)
samples = sampler.random(n=n_samples)
best_result = {'rmse': np.inf, 'parameters': {}}
for i, sample in enumerate(samples):
# Scale sample to parameter bounds
parameters = {}
for j, param_name in enumerate(param_names):
lower, upper = param_ranges[j]
parameters[param_name] = lower + sample[j] * (upper - lower)
# Simulate with these parameters
try:
simulated_ecg = self._simulate_with_parameters(parameters)
rmse = self.objective_function(parameters, simulated_ecg)
if rmse < best_result['rmse']:
best_result['rmse'] = rmse
best_result['parameters'] = parameters
if i % 100 == 0:
print(f"Sample {i}: RMSE = {rmse:.4f}%")
except Exception as e:
print(f"Error at sample {i}: {e}")
continue
self.optimal_parameters = best_result['parameters']
self.best_rmse = best_result['rmse']
return best_result
def _simulate_with_parameters(self, parameters: Dict[str, float]) -> Dict:
"""Simulate ECG with specified parameters"""
# Update model parameters
if 'san_location_alpha' in parameters:
# Modify stimulus site based on SAN location
san_alpha = parameters['san_location_alpha']
stimulus_idx = np.argmin(np.abs(self.model.uac['alpha'] - san_alpha))
else:
stimulus_idx = 0 # Default SAN location
# Solve activation
activation_time = self.model.solve_activation(stimulus_idx)
# Generate voltage time series
time_points = np.linspace(0, 0.15, 100)
v_m_series = []
for t in time_points:
v_m = self.model.compute_transmembrane_voltage(activation_time, t)
v_m_series.append(v_m)
v_m_series = np.array(v_m_series)
# Compute ECG
ecg_computer = LeadFieldECGComputation(
self.model.vertices, self.model.vertices,
self._get_electrode_positions()
)
simulated_ecg = ecg_computer.compute_ecg(v_m_series)
return simulated_ecg
def _get_electrode_positions(self) -> Dict:
"""Get standard 12-lead ECG electrode positions"""
# Simplified positions (relative to heart)
electrodes = {
'I': np.array([-1.0, 0.0, 0.0]),
'II': np.array([-0.5, -0.866, 0.0]),
'III': np.array([0.5, -0.866, 0.0]),
'aVR': np.array([-0.5, 0.289, 0.0]),
'aVL': np.array([-0.5, 0.289, 0.0]),
'aVF': np.array([0.0, -0.577, 0.0]),
'V1': np.array([0.0, 0.0, 0.3]),
'V2': np.array([0.0, 0.0, 0.4]),
'V3': np.array([0.0, 0.0, 0.5]),
'V4': np.array([0.0, 0.0, 0.6]),
'V5': np.array([0.0, 0.0, 0.7]),
'V6': np.array([0.0, 0.0, 0.8]),
}
return electrodes
# ============================================================================
# 9. INTEGRATION AND MAIN WORKFLOW
# ============================================================================
class AtrialsDigitalTwinWorkflow:
"""Complete automated end-to-end workflow for atrial digital twin generation"""
def __init__(self, ct_image: np.ndarray, voxel_size: float = 0.40,
target_mesh_resolution: float = 0.90):
self.ct_image = ct_image
self.voxel_size = voxel_size
self.target_resolution = target_mesh_resolution
self.results = {}
def run_complete_workflow(self) -> Dict:
"""Execute complete automated pipeline"""
print("=" * 60)
print("ATRIAL ELECTROPHYSIOLOGY DIGITAL TWIN GENERATION")
print("=" * 60)
# Stage 1: Segmentation
print("\n[1/8] Performing automated segmentation...")
segmentation = AutomatedSegmentation(self.ct_image, self.voxel_size)
cardiac_seg = segmentation.segment_by_intensity()
self.results['segmentation'] = cardiac_seg
# Stage 2: Mesh generation
print("[2/8] Generating volumetric meshes...")
mesh_gen = VolumetricMeshGenerator(cardiac_seg, self.voxel_size)
wall_volume = mesh_gen.extrude_atrial_walls(wall_thickness=3.0)
smoothed_wall = mesh_gen.smooth_surface(wall_volume)
la_vertices, la_faces = mesh_gen.generate_surface_mesh(smoothed_wall)
self.results['la_mesh'] = (la_vertices, la_faces)
# Stage 3: Landmark detection
print("[3/8] Detecting anatomical landmarks...")
landmark_detector = LandmarkDetection(la_vertices)
la_landmarks = landmark_detector.identify_la_landmarks()
self.results['la_landmarks'] = la_landmarks
# Stage 4: Universal Atrial Coordinates
print("[4/8] Computing Universal Atrial Coordinates...")
# Create dummy elements for UAC computation
n_verts = len(la_vertices)
dummy_elements = np.array([np.arange(4) + i for i in range(n_verts - 4)])
uac_computer = UniversalAtrialCoordinates(la_vertices, dummy_elements, 'LA')
la_uac = uac_computer.compute_uac()
self.results['la_uac'] = la_uac
# Stage 5: Fiber architecture
print("[5/8] Generating fiber architecture...")
fiber_gen = FiberArchitecture(la_vertices, la_uac)
fibers = fiber_gen.compute_fiber_orientation()
self.results['fibers'] = fibers
# Stage 6: EP Model setup
print("[6/8] Setting up electrophysiological model...")
ep_params = EPParameters(
cv_longitudinal=0.98,
cv_transverse=0.76,
conductivity_il=0.570,
conductivity_it=0.220,
conductivity_el=0.726,
conductivity_et=1.102,
surface_to_volume_ratio=1400.0
)
re_model = ReactionEikonalModel(la_vertices, dummy_elements, fibers, ep_params)
self.results['re_model'] = re_model
# Stage 7: Activation simulation
print("[7/8] Simulating atrial activation...")
stimulus_site = 0 # SAN location
activation = re_model.solve_activation(stimulus_site)
self.results['activation'] = activation
# Stage 8: ECG computation
print("[8/8] Computing 12-lead ECG...")
time_points = np.linspace(0, 0.15, 100)
v_m_series = []
for t in time_points:
v_m = re_model.compute_transmembrane_voltage(activation, t)
v_m_series.append(v_m)
v_m_series = np.array(v_m_series)
ecg_computer = LeadFieldECGComputation(la_vertices, la_vertices,
self._get_electrode_positions())
ecg = ecg_computer.compute_ecg(v_m_series)
self.results['ecg'] = ecg
print("\n" + "=" * 60)
print("WORKFLOW COMPLETED SUCCESSFULLY")
print("=" * 60)
return self.results
def _get_electrode_positions(self) -> Dict:
"""Standard 12-lead electrode positions"""
return {
'I': np.array([-1.0, 0.0, 0.0]),
'II': np.array([-0.5, -0.866, 0.0]),
'III': np.array([0.5, -0.866, 0.0]),
'aVR': np.array([-0.5, 0.289, 0.0]),
'aVL': np.array([-0.5, 0.289, 0.0]),
'aVF': np.array([0.0, -0.577, 0.0]),
'V1': np.array([0.0, 0.0, 0.3]),
'V2': np.array([0.0, 0.0, 0.4]),
'V3': np.array([0.0, 0.0, 0.5]),
'V4': np.array([0.0, 0.0, 0.6]),
'V5': np.array([0.0, 0.0, 0.7]),
'V6': np.array([0.0, 0.0, 0.8]),
}
# ============================================================================
# 10. EXAMPLE USAGE AND VISUALIZATION
# ============================================================================
if __name__ == "__main__":
# Create synthetic CT image
image_size = 256
ct_image = np.random.rand(image_size, image_size, image_size) * 0.5
# Create synthetic heart region
center = np.array([image_size // 2, image_size // 2, image_size // 2])
y, x, z = np.ogrid[:image_size, :image_size, :image_size]
# LA (ellipsoid)
la_mask = ((x - center[0])**2 / 25 + (y - center[1])**2 / 20 +
(z - center[2])**2 / 22) < 1
ct_image[la_mask] += 0.3
# RA (ellipsoid)
ra_mask = ((x - center[0] - 10)**2 / 20 + (y - center[1])**2 / 18 +
(z - center[2])**2 / 20) < 1
ct_image[ra_mask] += 0.25
print("Synthetic CT image created: shape =", ct_image.shape)
# Run complete workflow
workflow = AtrialsDigitalTwinWorkflow(ct_image, voxel_size=0.40,
target_mesh_resolution=0.90)
results = workflow.run_complete_workflow()
# Print results summary
print("\nRESULTS SUMMARY:")
print(f" - LA Mesh vertices: {results['la_mesh'][0].shape[0]}")
print(f" - LA Mesh faces: {results['la_mesh'][1].shape[0]}")
print(f" - LA UAC computed: {list(results['la_uac'].keys())}")
print(f" - Fiber orientations: {results['fibers'].shape}")
print(f" - Activation times (min/max): {results['activation'].min():.4f} / {results['activation'].max():.4f} s")
print(f" - ECG leads computed: {list(results['ecg'].keys())}")
# Plot sample results
plt.figure(figsize=(12, 8))
# Plot 1: Activation map
plt.subplot(2, 2, 1)
scatter = plt.scatter(results['la_mesh'][0][:, 0],
results['la_mesh'][0][:, 1],
c=results['activation'], cmap='viridis', s=1)
plt.colorbar(scatter, label='Activation Time (s)')
plt.title('Atrial Activation Map')
plt.xlabel('X (mm)')
plt.ylabel('Y (mm)')
# Plot 2: Sample ECG leads
plt.subplot(2, 2, 2)
for lead in ['I', 'II', 'III']:
if lead in results['ecg']:
plt.plot(results['ecg'][lead], label=lead, alpha=0.7)
plt.title('Sample ECG Leads (I, II, III)')
plt.xlabel('Time Sample')
plt.ylabel('Voltage (mV)')
plt.legend()
plt.grid(True, alpha=0.3)
# Plot 3: UAC Alpha coordinate
plt.subplot(2, 2, 3)
scatter = plt.scatter(results['la_mesh'][0][:, 0],
results['la_mesh'][0][:, 1],
c=results['la_uac']['alpha'], cmap='coolwarm', s=1)
plt.colorbar(scatter, label='α Coordinate')
plt.title('Universal Atrial Coordinate (α)')
plt.xlabel('X (mm)')
plt.ylabel('Y (mm)')
# Plot 4: Fiber orientation
plt.subplot(2, 2, 4)
sample_indices = np.random.choice(len(results['fibers']), 100, replace=False)
positions = results['la_mesh'][0][sample_indices]
directions = results['fibers'][sample_indices]
plt.quiver(positions[:, 0], positions[:, 1],
directions[:, 0], directions[:, 1], alpha=0.5)
plt.title('Fiber Orientation (Sample)')
plt.xlabel('X (mm)')
plt.ylabel('Y (mm)')
plt.tight_layout()
plt.savefig('atrial_digital_twin_results.png', dpi=150)
print("\nVisualization saved as 'atrial_digital_twin_results.png'")Related posts, You May like to read
- 7 Shocking Truths About Knowledge Distillation: The Good, The Bad, and The Breakthrough (SAKD)
- MOSEv2: The Game-Changing Video Object Segmentation Dataset for Real-World AI Applications
- MedDINOv3: Revolutionizing Medical Image Segmentation with Adaptable Vision Foundation Models
- SurgeNetXL: Revolutionizing Surgical Computer Vision with Self-Supervised Learning
- How AI is Learning to Think Before it Segments: Understanding Seg-Zero’s Reasoning-Driven Image Analysis
- SegTrans: The Breakthrough Framework That Makes AI Segmentation Models Vulnerable to Transfer Attacks
- Universal Text-Driven Medical Image Segmentation: How MedCLIP-SAMv2 Revolutionizes Diagnostic AI
- Towards Trustworthy Breast Tumor Segmentation in Ultrasound Using AI Uncertainty
- DVIS++: The Game-Changing Decoupled Framework Revolutionizing Universal Video Segmentation
- Radar Gait Recognition Using Swin Transformers: Beyond Video Surveillance
