[WIP] Reaction-Eikonal example

Project.toml

@@ -40,8 +40,8 @@ CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 GPUArrays = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
 
 [extensions]
-ThunderboltGPUArraysExt = "GPUArrays"
 CuThunderboltExt = ["CUDA", "GPUArrays"]
+ThunderboltGPUArraysExt = "GPUArrays"
 
 [compat]
 Aqua = "0.8"
@@ -72,6 +72,7 @@ julia = "1.10"
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+FastIterativeMethod = "51611a6c-7c63-4cd2-9a52-92f780132403"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 MatrixDepot = "b51810bb-c9f3-55da-ae3c-350fc1fbce05"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
@@ -84,5 +85,8 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
+[sources]
+FastIterativeMethod = {url = "https://git.noc.ruhr-uni-bochum.de/mohama17/eikonal.jl.git"}
+
 [targets]
-test = ["Aqua", "DelimitedFiles", "JET", "ModelingToolkit", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqTsit5", "Pkg", "SHA", "StaticArrays", "Tensors", "Test", "MatrixDepot"]
+test = ["Aqua", "DelimitedFiles", "JET", "ModelingToolkit", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqTsit5", "Pkg", "SHA", "StaticArrays", "Tensors", "Test", "MatrixDepot", "FastIterativeMethod"]

docs/Project.toml

@@ -1,7 +1,9 @@
 [deps]
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+Eikonal = "51611a6c-7c63-4cd2-9a52-92f780132403"
 Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992"
 FerriteGmsh = "4f95f4f8-b27c-4ae5-9a39-ea55e634e36b"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
@@ -13,7 +15,12 @@ OrdinaryDiffEqOperatorSplitting = "760fc936-9fa0-4281-9c1a-468957eedc89"
 OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 Thunderbolt = "909927c2-98d5-4a67-bba9-79f03a9ad49b"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
+
+[sources]
+Eikonal = {url = "https://git.noc.ruhr-uni-bochum.de/mohama17/eikonal.jl.git"}
+Thunderbolt = {path = ".."}

docs/make.jl

@@ -54,7 +54,7 @@ bibtex_plugin = CitationBibliography(
                 "EP02: Purkinje Network (TODO)" => "tutorials/ep02_purkinje.md",
                 "EP03: Defibrillation (TODO)" => "tutorials/ep03_bidomain.md",
                 "EP04: Monodomain ECG" => "tutorials/ep04_geselowitz-ecg.md",
-                "EP05: Eikonal Models (WIP)" => "tutorials/ep05_eikonal.md",
+                "EP05: Eikonal Models" => "tutorials/ep05_eikonal.md",
                 "EP06: Pacemakers (TODO)" => "tutorials/ep06_pacemaker.md",
             ],
         ],

docs/src/api-reference/models.md

@@ -101,6 +101,7 @@ MonodomainModel
 Thunderbolt.ParabolicParabolicBidomainModel
 Thunderbolt.ParabolicEllipticBidomainModel
 ReactionDiffusionSplit
+Thunderbolt.EikonalModel
 ```
 
 ```@docs

docs/src/assets/references.bib

@@ -360,6 +360,34 @@ @article{BakFalKolYan:2011:MSU
   abstract  = {This paper investigates the properties of smoothers in the context of algebraic multigrid (AMG) running on parallel computers with potentially millions of processors. The development of multigrid smoothers in this case is challenging, because some of the best relaxation schemes, such as the Gauss–Seidel (GS) algorithm, are inherently sequential. Based on the sharp two-grid multigrid theory from Falgout and Vassilevski (2004) and Falgout, Vassilevski, and Zikatanov (2005), we characterize the smoothing properties of a number of practical candidates for parallel smoothers, including several C-F, polynomial, and hybrid schemes. We show, in particular, that the popular hybrid GS algorithm has multigrid smoothing properties which are independent of the number of processors in many practical applications, provided that the problem size per processor is large enough. This is encouraging news for the scalability of AMG on ultraparallel computers. We also introduce the more robust ℓ₁ smoothers, which are always convergent and have already proven essential for the parallel solution of some electromagnetic problems.}
 }
 
+@article{NeicCamposPrassl:2017:ECE,
+title = {Efficient computation of electrograms and ECGs in human whole heart simulations using a reaction-eikonal model},
+journal = {Journal of Computational Physics},
+volume = {346},
+pages = {191-211},
+year = {2017},
+issn = {0021-9991},
+doi = {https://doi.org/10.1016/j.jcp.2017.06.020},
+url = {https://www.sciencedirect.com/science/article/pii/S0021999117304655},
+author = {Aurel Neic and Fernando O. Campos and Anton J. Prassl and Steven A. Niederer and Martin J. Bishop and Edward J. Vigmond and Gernot Plank},
+keywords = {Cardiac electrophysiology, Bidomain model, Eikonal model, Electrical activation and repolarization},
+abstract = {Anatomically accurate and biophysically detailed bidomain models of the human heart have proven a powerful tool for gaining quantitative insight into the links between electrical sources in the myocardium and the concomitant current flow in the surrounding medium as they represent their relationship mechanistically based on first principles. Such models are increasingly considered as a clinical research tool with the perspective of being used, ultimately, as a complementary diagnostic modality. An important prerequisite in many clinical modeling applications is the ability of models to faithfully replicate potential maps and electrograms recorded from a given patient. However, while the personalization of electrophysiology models based on the gold standard bidomain formulation is in principle feasible, the associated computational expenses are significant, rendering their use incompatible with clinical time frames. In this study we report on the development of a novel computationally efficient reaction-eikonal (R-E) model for modeling extracellular potential maps and electrograms. Using a biventricular human electrophysiology model, which incorporates a topologically realistic His–Purkinje system (HPS), we demonstrate by comparing against a high-resolution reaction–diffusion (R–D) bidomain model that the R-E model predicts extracellular potential fields, electrograms as well as ECGs at the body surface with high fidelity and offers vast computational savings greater than three orders of magnitude. Due to their efficiency R-E models are ideally suitable for forward simulations in clinical modeling studies which attempt to personalize electrophysiological model features.}
+}
+
+@article{Weingart:1977:AOIC,
+author = {Weingart, R},
+title = {The actions of ouabain on intercellular coupling and conduction velocity in mammalian ventricular muscle.},
+journal = {The Journal of Physiology},
+volume = {264},
+number = {2},
+pages = {341-365},
+doi = {https://doi.org/10.1113/jphysiol.1977.sp011672},
+url = {https://physoc.onlinelibrary.wiley.com/doi/abs/10.1113/jphysiol.1977.sp011672},
+eprint = {https://physoc.onlinelibrary.wiley.com/doi/pdf/10.1113/jphysiol.1977.sp011672},
+abstract = {1. The effects of ouabain on the electrical coupling between cells and the conduction velocity, theta, were studied in ventricular muscle preparations from calf and cow hearts using a silicon-oil-chamber. 2. After 90 min of exposure to 2 X 10(-6) M ouabain, an increase of the inside longitudinal resistance, Ri, from 420 omega cm to 1032 omega CM was observed. Assuming a constant myoplasmic resistivity this presumably reflects a reduced electrical coupling between myocardial cells. 3. Concomitantly, theta was decreased from 50-3 to 29-4 cm/sec. This change could be explained by the observed alterations in the maximal rate of rise of the action potential, (dV/dt)max, the amplitude of the action potential Vp, the membrane capacity Cf, and the sum, respectively, of the inside and outside longitudinal resistance per unit distance (ri + ro). Quantitatively, about 60\% of the decrease of theta could be accounted for by the experimentally determined increase of Ri. 4. Time course studies revealed a biphasic action of ouabain on Ri. An early dose-dependent drop in Ri, equivalent to an improvement of the intercellular coupling, was followed by a delayed massive increase in Ri, whose onset and magnitude were also concentration-dependent. 5. The delayed increase in Ri was associated with an increase of the diastolic tension. Toxic ouabain doses (2 X 10(-6) M) produced irreversible changes on both parameters, whereas thereapeutic doses (less than 5 X 10(-7) M) affected neither of them. Reversible effects on both parameters were observed at an intermediate drug concentration (10(-6) M). 6. The strong correlation between decoupling and contracture is consistent with the idea that the intracellular Ca concentration, [Ca]i, is involved in the control of the nexal conductance. This is supported by the finding that increasing the extracellular Ca concentration, [Ca]o, accelerated the ouabain-induced decoupling, whereas reducing [Ca]o retarded it. 7. If anything, the contracture slightly preceded the increase in Ri. From this it is concluded that the threshold [Ca]i for the electrical decoupling between cells must be somewhat larger than the threshold level for the tension activation. 8. The delayed increase in Ri is compatible with an inhibition of the Na pump which according to the Na-lag hypothesis predicts an increase of [Ca]i secondary to a Na-accumulation. The early drop in Ri can either be explained by a stimulation of the Na pump, or by a non-monotonic relationship between Ri and [Ca]i.},
+year = {1977}
+}
+
 @book{saad2003iterative,
   author    = {Saad, Yousef},
   title     = {Iterative Methods for Sparse Linear Systems},

docs/src/literate-tutorials/ep05_eikonal.jl

@@ -1,14 +1,183 @@
 # # [Electrophysiology Tutorial 5: Eikonal Models](@id ep-tutorial_eikonal)
 #
-# !!! todo
-#     Show activation timings.
+# ```@raw html
+# <img src="../eikonal-wave.gif" width="50%">
+# ```
 #
-# This tutorial shows how to solve Eikonal models and how to recover transemembrane potential fields.
+# This tutorial shows how to solve Eikonal models and how to recover transemembrane
+# potential fields and use them to calculate the ECG corresponding to the activation pattern.
 #
 # !!! todo
 #     Provide context.
 #
+# # Introduction
+# Depolarization can be modeled spatially as a wave propagating through cardiac tissue.
+# The arrival times of this wave is the solution of the Eikonal equation in the Reaction-Eikonal
+# split [NeicCamposPrassl:2017:ECE](@citet):
+# ```math
+# \begin{aligned}
+#   C_{\textrm{m}} \partial_t \varphi &=  I_{\textrm{foot}}(t_{\textrm{arrival}}, t) - I_{\textrm{ion}}(\varphi, \boldsymbol{s}, t) & \textrm{pointwise in} \: \Omega \, , \\
+#   \partial_t \boldsymbol{s} &= \mathbf{g}(\varphi, \boldsymbol{s}) & \textrm{pointwise in}  \: \Omega \, , \\
+#   I_{\textrm{foot}}(t_{\textrm{arrival}}, t) &= -\frac{A}{\tau}\cdot e^{\frac{t - t_{\textrm{arrival}}}{\tau}} & \textrm{pointwise in}  \: \Omega \;\,
+# \end{aligned}
+# ```
+# Where $C_{\textrm{m}}$ is the membrane capacitance, $\varphi$ is the transmembrane
+# potential, $I_{\textrm{foot}}$ is the stimulating current, $I_{\textrm{ion}}$ is
+# the ionic cureent resulting from the cell model. $\boldsymbol{s}$ is the state vector
+# for the cell model, and $\mathbf{g}$ is the function evolving the cell model states in time.
+#
+# Wave arrival times are obtained by solving the Eikonal equation:
+# ```math
+# \sqrt{\nabla t_{\textrm{arrival}}^T \boldsymbol{V} \nabla t_{\textrm{arrival}}} = 1 \qquad \textrm{in}  \: \Omega
+# ```
+#
+# Where $\boldsymbol{V}$ the conduction velocity, which was shown to be proportional to
+# the square root of the diffusivities of the cardiac tissue [Weingart:1977:AOIC](@citet) (That's the earliest mention I
+# could find but I remember there was an older one in the 50s).
+#
+#
 # ## Commented Program
+using Thunderbolt, LinearAlgebra, StaticArrays, DifferentialEquations, FastIterativeMethod
+
+# !!! todo
+#     The initializer API is not yet finished and hence we deconstruct stuff here manually.
+#     Please note that this method is quite fragile w.r.t. to many changes you can make in the code below.
+function steady_state_initializer!(u₀, f::Thunderbolt.EikonalCoupledODEFunction)
+    ## TODO cleaner implementation. We need to extract this from the types or via dispatch.
+    odefun = f.ode_function.prob_func
+    ionic_model = odefun.cellmodel
+
+    φ₀ = @view u₀[1:solution_size(f.eikonal_function)]
+    ## TODO extraction these via utility functions
+    s₀flat = @view u₀[(solution_size(f.eikonal_function) + 1):end]
+    ## Should not be reshape but some array of arrays fun
+    s₀ = reshape(
+        s₀flat, (solution_size(f.eikonal_function), Thunderbolt.num_states(ionic_model) - 1)
+    )
+    default_values = Thunderbolt.default_initial_state(ionic_model)
+
+    φ₀ .= default_values[1]
+    for i in 1:(Thunderbolt.num_states(ionic_model) - 1)
+        s₀[:, i] .= default_values[i + 1]
+    end
+    return
+end
+
+# We start by generating a tetrahedral mesh, as the fast iterative method solver
+# accepts tetrahedral meshes only.
+heart_mesh = generate_ideal_lv_mesh(
+    12, 2, 5;
+    inner_radius = 0.7,
+    outer_radius = 1.0,
+    longitudinal_upper = 0.0,
+    apex_inner = 0.7,
+    apex_outer = 1.0
+) |> Thunderbolt.hexahedralize |> Thunderbolt.tetrahedralize;
+
+# !!! tip
+#     We can also load realistic geometries with external formats. For this simply use either FerriteGmsh.jl
+#     or one of the loader functions stated in the [mesh API](@ref mesh-utility-api).
+
+# We also generate a coordinate system to be used for checking points for initial activation.
+cs = compute_lv_coordinate_system(heart_mesh)
+
+# For this tutorial we define a uniform endocardial activation
+activation_protocol = Thunderbolt.UniformEndocardialEikonalActivationProtocol()
+
+# For our toy problem we use a very simple microstructure.
+microstructure = OrthotropicMicrostructureModel(
+    ConstantCoefficient((Vec(0.0, 0.0, 1.0))),
+    ConstantCoefficient((Vec(0.0, 1.0, 0.0))),
+    ConstantCoefficient((Vec(1.0, 0.0, 0.0)))
+)
+
+# With the microstructure we setup the diffusion tensor field in spectral form.
+# !!! todo
+#     citation
+κ₁ = 0.17 * 0.62 / (0.17 + 0.62)
+κᵣ = 0.019 * 0.24 / (0.019 + 0.24)
+diffusion_tensor_field = SpectralTensorCoefficient(
+    microstructure,
+    # ConstantCoefficient(SVector(κ₁, κᵣ, κᵣ))
+    ConstantCoefficient(SVector(κ₁, κ₁, κ₁))
+)
+# Now we choose our *inner* cell model, used to compute the ionic currents, since
+# the eikonal-based foot current acts as a wrapper around ionic cell models.
+cellmodel = Thunderbolt.PCG2019()
+
+# We define the monodomain parameters to be used in the split.
+heart_model = MonodomainModel(
+    ConstantCoefficient(1.0),
+    ConstantCoefficient(1.0),
+    diffusion_tensor_field,
+    NoStimulationProtocol(),
+    cellmodel,
+    :φₘ, :s
+)
+
+# We semidiscretize the using a reaction-diffusion split and provide the desired activation
+# protocol. Here, we activate the full indocardium defined as nodes with less than 0.05
+# transmural coordinate in the provided coordinate system by default.
+heart_odeform = semidiscretize(
+    ReactionEikonalSplit(heart_model, cs),
+    (
+        FiniteElementDiscretization(Dict(:φₘ => LagrangeCollection{1}())),
+        Thunderbolt.SimplicialEikonalDiscretization(;
+            activation_protocol
+        ),
+    ),
+    heart_mesh
+)
+
+# We allocate the solution vector for the reaction part. The activation times vector is allocated internally.
+u₀ = zeros(Float64, length(heart_odeform.ode_function.prob.u0) * length(heart_mesh.grid.nodes))
+
+# Apply the initial conditions.
+steady_state_initializer!(u₀, heart_odeform)
+
+# Here we solve for the wave arrival time using the package *Eikonal.jl* as it's
+# required to be solved for once before timestepping.
+FastIterativeMethod.solve!(heart_odeform, heart_mesh, diffusion_tensor_field)
+
+# !!! danger
+#     The fast iterative method solver in `Eikonal.jl` uses `@inbounds` a lot
+#     in its internals. It was observed that forcing bounds checks results in an
+#     enormous amount of dynamic allocations that can easily result in running out of memory.
+
+nodal_timings = heart_odeform.ode_function.prob_func.activation_timings
+VTKGridFile(("nein-activation-map-debug.vtu"), heart_mesh.grid) do vtk              # hide
+    vtk.vtk["timings"] = nodal_timings    # hide
+end                                                                                 # hide
+
+dt₀ = 0.01
+dtvis = 0.5
+Tₘₐₓ = 50.0
+Tₘₐₓ = dtvis # hide
+tspan = (0.0, Tₘₐₓ)
+
+
+sim = solve(heart_odeform.ode_function, Rodas5P(), saveat = collect(0.0:dtvis:Tₘₐₓ), trajectories = length(heart_odeform.ode_function.prob_func.activation_timings))
+dh = cs.dh |> deepcopy
+Thunderbolt.reorder_nodal!(dh)
+
+io = ParaViewWriter("ep05_eikonal")
+
+φₘfield = copy(nodal_timings)           # hide
+for t in 0.0:dtvis:Tₘₐₓ                                                          # hide
+    for i in 1:length(φₘfield)                                                  # hide
+        φₘfield[i] = sim.u[i](t)[1]                                             # hide
+    end                                                                         # hide
+    store_timestep!(io, t, dh.grid) do file                                     # hide
+        Thunderbolt.store_timestep_field!(file, t, dh, φₘfield, :coordinates)   # hide
+    end                                                                         # hide
+end                                                                             # hide
+
+
+## test the result                                                      #src
+using Test                                                              #src
+@test maximum(nodal_timings) ≈ 0.3 / κ₁                                   #src
+@test count(≈(0.3 / κ₁), nodal_timings) == (5 + 1) * 12 + 1                     #src
+@test count(≈(0.3 / κ₁ * cosd(90 / 12)), nodal_timings) == (5 + 1) * 12           #src
 
 
 #md # ## References

docs/src/tutorials/index.md

@@ -93,7 +93,7 @@ In this tutorial you will learn how to:
 * transfer coefficients and solutions between overlapping domains
 * compute the ECG form a monodomain model as a postprocessing step
 ---
-#### [EP Tutorial 05: Reaction-Eikonal ECG](@ref ep-tutorial_eikonal) (TODO)
+#### [EP Tutorial 05: Reaction-Eikonal ECG](@ref ep-tutorial_eikonal)
 In this tutorial you will learn how to:
 * perform ECG simulations with a simplified activation dynamics model
 ---

src/Thunderbolt.jl

@@ -88,6 +88,9 @@ include("modeling/solid_mechanics.jl")
 include("modeling/fluid_mechanics.jl")
 
 include("modeling/multiphysics.jl")
+include("modeling/core/eikonal.jl")
+
+include("discretization/eikonal.jl")
 
 include("modeling/functions.jl")
 include("modeling/problems.jl")
@@ -198,6 +201,8 @@ export
     TransmembraneStimulationProtocol,
     AnalyticalTransmembraneStimulationProtocol,
     ReactionDiffusionSplit,
+    ReactionEikonalSplit,
+    EikonalCoupledODEFunction,
     # Circuit
     RSAFDQ2022LumpedCicuitModel,
     MTKLumpedCicuitModel,