Hi, I'm trying to simulate ventricular electrophysiology on a 3D mesh. My 2D tissue simulations work perfectly with realistic wave propagation, but the same setup on 3D ventricular meshes produces no activation (all LATs = -1.0) with clinically-relevant stimulus strengths, or simultaneous activation (no propagation) with extremely high strengths.
My setup is as follows:
- Mesh: Extracted left and right ventricles (tags 1 and 2) from the Strocchi et al. four-chamber mesh (Zenodo 3890034), then resampled to ~2mm resolution (114,261 nodes, 588,654 elements)
- Model: Ten Tusscher-Panfilov with modified GKs conductance
- Solver: Monodomain (bidomain shows same behavior)
- Electrode: 82 nodes from clinical RV endocardial pacing site, remapped after resampling using nearest-neighbor matching
- Conductivities: g_il=0.174, g_it=0.019, g_in=0.019 S/m (anisotropic, from working 2D example)
What Worked:
2D tissue simulations (10,000 nodes, from openCARP example 02_EP_tissue/08_lats):
- Pulse strength: 100 μA/cm²
- Result: Realistic propagation with ACTs spreading from 10ms to 21ms
- Full parameters in working case available on request
What Didn't Work:
3D ventricular mesh with identical ionic model and conductivity parameters:
| Pulse Strength | Result |
|---|
| 100 (clinical) | All ACTs = -1.0 |
| 500 | All ACTs = -1.0 |
| 50,000 | All ACTs identical within 0.002ms (simultaneous, no propagation) |
Specifically, I took these steps to get the 3D meshes that I worked with:
- Downloaded Strocchi mesh 01.tar.gz from Zenodo 3890034
- Converted from EnSight Gold (.case) to openCARP using: meshtool convert -imsh 01-350um.vtk -ifmt vtk_bin -omsh 01-350um -ofmt carp_txt
- Extracted ventricles: meshtool extract mesh -tags 1,2
- Resampled to 2mm: meshtool resample mesh -min 1.5 -max 2.5
- Remapped electrode using scipy.spatial.cKDTree nearest-neighbor matching
Also, even without extracting ventricles, I subsampled to a smaller box within the original (converted) mesh (i.e. 01-350um) for faster testing, and was getting similar issues, so I don't think the problem was with extracting ventricles or resampling to 2mm. This is how I got the smaller box:
meshtool extract volume \
-msh /workspace/heart_3d/01_convert/01-350um \
-ifmt carp_txt \
-submsh /workspace/heart_3d/01_tiny_full/tiny_block \
-ofmt carp_txt \
-coord 40,2,100,2,360,50
Playing around with the size of the box and the size of the stimulus box (defined by stim[0].elec.p0[0], p0[1], p0[2], p1[0], etc), I noticed that it seems to relate to stimulus strength: the larger the box i choose, the higher the pulse strength needs to be for activation to happen (for example, with a 3x3x9mm box, which comprised 1,448 nodes, pulse strength of 50 was sufficient to produce activations – although, again, they were all nearly identical and did not propagate spatially; whereas with a 10x10x10 box of 24,000 nodes, i needed a pulse strength of 50000 or so to get any activations in the lats file, which were again nearly identical).
So across all the things I've tried with the 3d mesh, even when activation occured, all nodes activate simultaneously (no spatial propagation)
I also performed some diagnostics:
Fiber orientations: 430,555 unique vectors (appropriate variety)
Mesh connectivity: Single connected component confirmed via meshtool
Electrode remapping: Verified nearest-neighbor remapping from original to resampled mesh
Stimulus firing: Stimulus_0.trc file shows current pulse applied correctly
Domain specification: Tried both "intra" and "extra" in .vtx file
Physics: Tested both monodomain and bidomain
Isotropy: Tested isotropic conductivity (g_il=g_it=g_in=0.5) - no improvement
Is there a recommended workflow for going from the Strocchi .case/.geo format to working openCARP simulations? Or alternatively, is there a known issue with this dataset that requires special handling, and is there a better dataset that I should be using for ventricular 3D heart simulations?
Here is an example of code that I used for the runs:
mpiexec -n 24 openCARP \
-meshname ventricles_2mm \
-ellip_use_pt 0 \
-parab_use_pt 0 \
-parab_options_file /usr/local/lib/python3.10/distpackages/carputils/resources/petsc_options/ilu_cg_opts \
-ellip_options_file /usr/local/lib/python3.10/distpackages/carputils/resources/petsc_options/gamg_cg_opts \
-num_phys_regions 2 \
-phys_region[0].ptype 0 \
-phys_region[0].num_IDs 2 \
-phys_region[0].ID[0] 1 \
-phys_region[0].ID[1] 2 \
-phys_region[1].ptype 1 \
-phys_region[1].num_IDs 2 \
-phys_region[1].ID[0] 1 \
-phys_region[1].ID[1] 2 \
-num_gregions 2 \
-gregion[0].num_IDs 1 \
-gregion[0].ID[0] 1 \
-gregion[0].g_il 0.174 \
-gregion[0].g_it 0.019 \
-gregion[0].g_in 0.019 \
-gregion[0].g_el 0.625 \
-gregion[0].g_et 0.236 \
-gregion[0].g_en 0.236 \
-gregion[1].num_IDs 1 \
-gregion[1].ID[0] 2 \
-gregion[1].g_il 0.174 \
-gregion[1].g_it 0.019 \
-gregion[1].g_in 0.019 \
-gregion[1].g_el 0.625 \
-gregion[1].g_et 0.236 \
-gregion[1].g_en 0.236 \
-num_imp_regions 2 \
-imp_region[0].num_IDs 1 \
-imp_region[0].ID[0] 1 \
-imp_region[0].im tenTusscherPanfilov \
-imp_region[0].im_param "GKs*1.0" \
-imp_region[1].num_IDs 1 \
-imp_region[1].ID[0] 2 \
-imp_region[1].im tenTusscherPanfilov \
-imp_region[1].im_param "GKs*1.0" \
-num_stim 1 \
-stim[0].crct.type 0 \
-stim[0].pulse.strength 100.0 \
-stim[0].ptcl.duration 2.0 \
-stim[0].ptcl.npls 1 \
-stim[0].elec.vtx_file electrode_endo_rv \
-num_LATs 2 \
-lats[0].ID ACTs \
-lats[0].all 0 \
-lats[0].mode 0 \
-lats[0].threshold -10 \
-lats[1].ID REPs \
-lats[1].all 0 \
-lats[1].mode 1 \
-lats[1].threshold -70 \
-dt 0.025 \
-tend 400 \
-bidomain 0 \
-simID "3d_heart"
This particular run resulted in all -1.0's in the activation file. The electrode file looked like this:
82
extra
8370
9953
12497
... (etc)
Let me know if this makes sense and if there is more info I can provide. Maybe there is something obviously wrong with my command – I am new to openCARP. Thank you!
