Introduction

Human heart is an intrinsic multiphysics system. It’s normal function requires synchronization and coupling of cardiac electrophysiology, tissue (myocardium) mechanics, and hemodynamics. Among the three physics, cardiac mechanics plays a central role, as it serves as the bridge between electrophysiology and blood flow. In a physiological setting, the myocardial cells (myocytes) in the heart tissue initiate and propagate an electrical impulse (action potential) that travels through the heart and drives the contraction of the cardiac muscle, which, in turn, propels the blood to sustain whole body circulation.

The SimVascular SimCardio software suite includes a machine learning-based tool to automatically segment heart chambers from medical data and create patient-specific cardiac models, a finite element solver for performing multiphysics cardiac simulations, and python-based plugins to interactively create cardiac models, Purkinje network and set solver parameters from within the main SimVascular GUI application. While some of the interactive plugins are under development, the software allows flexibility to use scripting and perform simulations on high performance clusters.

Cardiac Segmentation and Geometric Modeling

The SimVascular Cardiac Image Segmentation Tool uses pre-trained deep neural network models to automaticall generate segmentations for major cardiac structures, the four heart chambers, aorta and pulmonary arteries. The automatic cardiac segmentation tool utilized an ensemble of two-dimensional (2D) convolutional neural networks (CNNs) for automatic segmentation of cardiac structures from three-dimensional (3D) patient images and demonstrated state-of-the-art performance than prior approaches when evaluated on a benchmark dataset containing both magnetic resonance (MR) and computed tomography(CT) scans. This tool can be used from both the SimVascular’s Python plugin in the GUI or from the command line version of the SimVascular Python.

Multiphysics Finite Element Solver

svFSI is a new multi-physics finite element solver designed for computational modeling of the whole heart dynamics. As the next generation underlining solver for the SimVascular software, svFSI is capable of simulating myocardial electrical activity using state-of-the-art cellular activation models, employ nonlinear hyperelastic material models to simulate tissue mechanics and activation coupling strategies, and perform large deformation fluid-structure interaction (FSI) to simulate the dynamics between the heart tissue and blood flow. svFSI could also be used to perform blood flow simulations by imposing the ventricular wall motion extracted from medical image data.

svFSI is compatible with several mesh formats and can read a wide range of 2D and 3D element types. svFSI could be coupled to Trilinos linear solver library providing a wide choice of linear solvers and preconditioners for various applications. The methods employed are widely popular among the scientific community and continue to evolve through research and development. The solver is parallelized using MPI for inter-process communication and was demonstrated to scale efficiently on large supercomputing clusters. ParMETIS is used for partitioning the computational domains. Simulation results are outputted into VTK format that can be easily visualized in the free software Paraview.

The source code of svFSI is available on GitHub as part of the SimVascular project and the precompiled installers can be downloaded from SimTK website. A variety of examples are also available on GitHub for the users to get started quickly on using the solver. We recommend users to navigate through all the examples under a particular physics to educate themselves with the variety of options that the solver could accommodate. We also recommend users to use the existing input files as templates and make suitable modifications that meet their needs.

Cardiac Geometric Modeling Tool

The SimVascular Cardiac Geometric Modeling Tool is used to extract major cardiac structures from time-resolved imaging data using automatic machine learning techniques. Volumetric meshes created from segmentation geometry are registered with image data to extract heart wall motion that can be later used in multi-physics simulations of the heart using the svFSI solver.

The Cardiac Geometric Modeling Tool is implemented as the sv_auto_lv_modeling Python package included with the SimVascualar application.

Automated Image Segmentation

The SimVascular Cardiac Geometric Modeling Tool uses pre-trained deep neural network models to automatically generate segmentations of major cardiac structures: the four heart chambers, aorta and pulmonary arteries. The segmentation tool uses an ensemble of two-dimensional (2D) convolutional neural networks (CNNs) for automatic segmentation of cardiac structures from three-dimensional (3D) patient images. It has been shown to have superior performance compared to prior approaches when evaluated on a benchmark dataset containing both magnetic resonance (MR) and computed tomography (CT) image data [1]. The tool can be used from both the SimVascular GUI Python console and from the command line within the SimVascular Python environment.

Input Requirements

The preferred input format of the image volumes is .nii.gz, .nii or .vti files. The segmentation method requires an image to have an identity orientation matrix.

The VTK version used by SimVascular does not writte a .vti file with an orientation matrix. Image .vti files written by SimVascular must therefore be reoriented to have an identity orientation matrix when used for segmentation.


The directory containing the input image data should be organized as follows

image_dir
     |__ patient_id (optional)
         |__ image_volume0.vti
         |__ image_volume1.vti
         |__ image_volume2.vti
         |__ ...

Pre-Trained Models

We used the image and ground truth data provided by MMWHS to train our models. Our segmentation models were trained simultaneously on CT and MR data. The trained weight data can be downloade from here.

Segmenting Image Data

Segmenting using the SimVascular Python shell

3D CT or MR image volumes can be segmented within the SimVascular Python Shell using the prediction.py Python script. For example


    data_path=/path/to/data
    sv_python_dir=/usr/local/bin
    script_dir=SimVascular/Python/site-packages/sv_auto_lv_modeling

    patient_id=WS01
    image_dir=$data_path/01-Images
    output_dir=$data_path/02-Segmnts
    weight_dir=$data_path/Weights

    ${sv_python_dir}/simvascular --python -- $script_dir/segmentation/prediction.py \
        --pid $patient_id \
        --image $image_dir \
        --output $output_dir \
        --model $weight_dir \
        --view  0 1 2 \ # 0 for axial, 1 for coronal, 2 for sagital
        --modality ct # ct or mr

Segmenting using the SimVascular GUI Python Console

Image segmentation can also be performed from within the SimVascular GUI Python console . Use the Python console Text Editor mode to enter the following Python commands


    from auto_lv.auto_lv import Segmentation
    data_path='/path/to/data'
    seg = Segmentation()
    seg.set_modality('ct')
    seg.set_patient_id ('WS01')
    seg.set_image_directory (data_path+'/01-Images')
    seg.set_output_directory (data_path+'/02-Segmnts')
    seg.set_model_directory ([data_path+'/Weights'])
    seg.set_view ([2])
    seg.generate_segmentation()

and execute them.

Visualizing Segmentation Results

Segmentation results can be visualized using SimVascular by first creating a SimVascular project. The segmentation results .vti file is then added as an image by right-clicking Images under the SV Data Manager and selecting the Add/Replace Image menu item. The following two pictures show an example of CT image and the segmentation results.

CT image data.
Segmentation results of the CT image data.

Segmentation results can be visualized using ParaView.

Mesh Generation

The Cardiac Geometric Modeling Tool can be used to generate left ventricle (LV) meshes from segmentation results. These meshes can then be used CFD modeling of LV flow combined with patient medical imaging data. The LV Mesh Generation tool is used to automatically generate meshes for simulation. It takes as input a segmentation and outputs reconstructed LV surface meshes.

Python scripts are provided for extracting a series of surfaces representing heart wall motion by registering surface meshes to time-resolved imaging data using the SimpleElastix medical image registration software. These time-series surface meshes can also be interpolated temporally to obtain the resolution needed for simulations.

Constructing LV Surface Meshes

Using SimVascular Python Shell

The SimVascular Python Shell can be used to run the surface_main.py Python surface generation script using the following commands


  # Path to SimVascular executable
  data_path=/path/to/data
  sv_python_dir=/usr/local/bin
  script_dir=SimVascular/Python/site-packages/sv_auto_lv_modeling
  # Path to the segmentation results
  p_id=WS01
  input_dir=$data_path/02-Segmnts/$p_id
  # Path to the outputed surface meshes
  output_dir=$data_path/03-Surfaces/$p_id
  model_script=$script_dir/modeling/surface_main.py
  # Construct LV surface meshes with tagged boundary faces
  ${sv_python_dir}/simvascular --python \
      -- ${model_script} \
      --input_dir ${input_dir} \
      --output_dir ${output_dir} \
      --edge_size 3.5 # maximum edge size for the mesh

Using SimVascular Python Console

The SimVascular GUI Python Console can be used for surface generation. Use the Python console Text Editor mode to enter the following Python commands

from sv_auto_lv_modeling.auto_lv import Modeling
    data_path='/path/to/data'
    surf = Modeling()
    surf.set_segmentation_directory(data_path+'/02-Segmnts/WS01')
    surf.set_output_directory (data_path+'/03-Surfaces/WS01')
    surf.set_max_edge_size (3.5)
    surf.generate_lv_modes ()
Automatically created LV meshes at diastole and systole.

Volumetric Meshing

Using SimVascular Python Shell

The volumetric meshing script volume_mesh_main.py is used to generate a finite element mesh.


    # Path to SimVascular exectuable
    data_path=/path/to/data
    sv_python_dir=/usr/local/bin
    script_dir=SimVascular/Python/site-packages/sv_auto_lv_modeling
    p_id=WS01

    # Path to the surface meshes
    input_dir=$data_path/04-SurfReg/$p_id

    # Path to the outputed volume meshes
    output_dir=$data_path/05-VolMesh/$p_id
    volume_mesh_script=$script_dir/modeling/volume_mesh_main.py

    # Volumetric Meshing using SimVascular
    ${sv_python_dir}/simvascular --python \
        -- ${volume_mesh_script} \
        --input_dir $input_dir \
        --output_dir $output_dir \
        --phase 0 \ # the phase id in $input_dir to generate a volumetric mesh
        --edge_size 3.5

Using the SimVascular Python Console

Within the Python plugin, we can use the Text Editor mode and enter the following lines to create a Python script.


    from sv_auto_lv_modeling.auto_lv import VolumeMesh
    data_path='/path/to/data'
    vol = VolumeMesh()
    vol.set_output_directory (data_path+'/05-VolMesh/WS01')
    vol.set_max_edge_size (3.5)
    vol.set_surface_model_filename (data_path+'/04-SurfReg/WS01/WS01_0.vti.vtp')
    vol.generate_volume_mesh()

Mesh Registration

We can simulate the LV flow over time by tracking the deformation of the heart from time-resolved imaging and impose this motion to the fluid domains inside the heart, which leads to a deforming-domain CFD problem. To be able to track the deformation of the generated LV mesh over time, we need to building point-corresponded LV meshes from segmentations at all time frames. We generate a surface mesh at one time frame and propagated to the others using the displacement field obtained from registering the corresponding segmentations.

The elastix_main.py script is used to perform the registration. It uses uses the SimpleElastix software for registration.

For code compatiblty you will need to use SimpleElastix commit 8244e0001f4137514b0f545f1e846910b3dd7769.



    # Use SimpleElastix to register surface meshes
    data_path=/path/to/data
    sv_python_dir=/usr/local/bin
    script_dir=SimVascular/Python/site-packages/sv_auto_lv_modeling

    # Path to the ct/mr images or segmentation results
    p_id=WS01
    image_dir=$data_path/01-Images/$p_id
    mask_dir=$data_path/02-Segmnts/$p_id

    # Path to the unregistered surface mesh
    surface_dir=$data_path/03-Surfaces/$p_id

    # Path to the registered surface meshes
    output_dir=$data_path/04-SurfReg/$p_id

    # Phase ID of the surface mesh used as the registration source
    start_phase=0

    # Registration with SimpleElastix
    python $script_dir/modeling/elastix_main.py \
        --image_dir $mask_dir \
        --mask_dir $mask_dir \
        --output_dir $output_dir \
        --start_phase $start_phase \
        --surface_dir $surface_dir \
        --image_file_extension vti \
        --edge_size 3.5

Mesh Motion

Once meshes have been registered the displacement of each mesh vertex over the whole cardiac cycle can be computed. The series of registered meshes are interpolated using cubic splines to obtain mesh displacements with a temporal resolution suitable for simulations.

The interpolation.py script is used to interpolate meshes and compute the mesh motion. It writes out a .dat file for each boundary face that can used in svFSI to set up the displacement boundary conditions


    # Generate motion.dat File for svFSI
    data_path=/path/to/data
    sv_python_dir=/usr/local/bin
    script_dir=SimVascular/Python/site-packages/sv_auto_lv_modeling

    # Phase ID should be the same as the one used in volume_mesh.sh
    phase_id=0
    p_id=WS01

    # Path to the registered surface meshes
    input_dir=$data_path/04-SurfReg/$p_id

    # Path to the outputted volumetric meshes
    output_dir=$data_path/05-VolMesh/$p_id

    # Number of interpolations between adjacent phases
    num=99

    # Number of cardiac cycles
    cyc=1

    # Cycle duration in seconds
    period=1.25

    # Write boundary conditions for FSI simulations
    python $script_dir/modeling/svfsi/interpolation.py \
        --input_dir $input_dir \
        --output_dir $output_dir \
        --num_interpolation $num \
        --num_cycle $cyc \
        --duration $period \
        --phase $phase_id





References

[1] Kong, F., and Shadden, S. C. (2020). Automating Model Generation for Image-based Cardiac Flow Simulation. ASME. J Biomech Eng.






Cardiac Electrophysiology Simulation Guide

Introduction

The electrical activity in a heart is initiated at the sinoatrial node, which is located at the right atrium of the heart. A small cluster of pacemaker cells periodically emits electrical signals into the heart’s conduction system. These signals first travel at a high speed through a network of special cells that allows the signals to reach the whole heart rapidly. Finally, they travel into the ventricular muscles, and cause the entire heart to contract almost simultaneously. The electrical signals travel between cells through a process called depolarization, during which voltage-sensitive protein channels on the cell membrane open to allow positively charged ions to move in and out of the cell, leading to ionic currents.

Mathematical Model of Cardiac Electrophysiology

The propagation of electrical signal in the heart is governed by a reaction-diffusion like equation,

$$\frac{\partial V}{\partial t} + \frac{I_{ion} - I_{app}(t)}{C_m} = \nabla \cdot\left( \mathbf{D}\nabla V \right) $$ $$ \mathrm{in~} \Omega^E\times(0,T] \nonumber$$ $$(\mathbf{D}\nabla V) \cdot \mathbf{N}=0 \mathrm{~on~} \partial\Omega^E\times(0,T]$$

$V$ is the local transmembrane potential. $C_m$ is the local membrane capacitance per unit area. $I_{ion}$ and $I_{app}$ are the ionic current flux (current per unit area) and applied external current flux, respectively. Here, $\mathbf{D}$ dictates the propagation velocity of the electrical signal and has the similar physical meaning as the diffusivity. It is calculated as

$$ \mathbf{D} =\frac{\sigma}{\chi C_m} $$

$\sigma$ is the conductivity tensor and has the unit $S/m$. $\chi$ is the ratio of membrane area over tissue volume. Then, $\mathbf{D}$ is a tensor with the dimension $Length^2/Time$.

The above equation is the mono-domain description of the cardiac electrophysiology, i.e. we don’t solve the intra- and extra-cellular electrical signal propagation separately. The mono-domain and multi-domain conductivities are connected through the following relation,

$$ \sigma = \frac{\sigma_i\sigma_e}{\sigma_i + \sigma_e} $$

where $\sigma_i$ and $\sigma_e$ are the intra- and extra-cellular conductivity tensors [1]. It is commonly assumed that the conductivity is transversely isotropic,

$$ \mathbf{\sigma} = \sigma_f \mathbf{f}\otimes \mathbf{f} + \sigma_s (\mathbf{I}-\mathbf{f}\otimes \mathbf{f}) $$

where $\sigma_f$ and $\sigma_s$ are the conductivities along the fiber direction and in the transverse plane. $\mathbf{f}$ is the fiber direction vector.

In svFSI, we directly specify $\mathbf{D}$ in the input file. We choose a slightly different form to enforce the transverse isotropy of the conductivity tensor,

$$ \mathbf{D} = D_{iso}\mathbf{I} + \sum_{n=1}^{nsd}D_{ani,n}\mathbf{fN}_n\otimes\mathbf{fN}_n $$

Here, $nsd$ is the dimension, and $\mathbf{fN_n}$ is the local orthonormal coordinate system built by fiber direction and sheet direction. To connect with the previous expression, we have $D_f=D_{iso}+D_{ani}$ and $D_s=D_{iso}$.

Cellular Activation Models

Depending on how the depolarization and repolarization within a single cardiac myocyte is described, the electrophysiology models can be roughly divided into two categories: biophysics-based ionic models (such as the ten Tusscher-Panfilov (TTP) model[2][3]), and phenomenological models (such as the Aliev-Panfilov (AP), Fitzhugh-Nagumo (FN) models[4]).

Biophysics-based Ionic Models

Between cardiac myocytes, the propagation of the electrical signal is enabled by the transmembrane motion of different ions, such as $K^+$, $Na^+$ and $Ca^{2+}$ during depolarization and repolarization of the myocytes. Cellular biophysics-based activation models are designed to capture these ionic movements. One of the popular biophysics-based ionic models, the TTP model [3] is implemented in svFSI. Here, we will briefly illustrate the ionic currents involved in the TTP model.

In the TTP model, the ionic current is expressed as

$$ I_{ion} = I_{Na} + I_{Kl} + I_{to} + I_{Kr} + I_{Ks} + I_{CaL} + I_{NaCa} + I_{NaK} + I_{pCa} + I_{pK} + I_{bCa} + I_{bNa} $$

The governing equations for each current can be found in [3].

Since, intracellular calcium concentration is the driving factor behind excitation-contraction coupling, we focus on the calcium dynamics here. In the TTP model,the following calcium concentrations are included:

  • $Ca_{itotal}$: total (free+buffered) cytoplasmic $Ca^{2+}$ concentration.
  • $Ca_{SRtotal}$: total SR $Ca^{2+}$ concentration.
  • $Ca_{SStotal}$: total dyadic space $Ca^{2+}$ concentration.
  • $Ca_i$: free cytoplasmic $Ca^{2+}$ concentration.
  • $Ca_{SR}$: free SR $Ca^{2+}$ concentration.
  • $Ca_{SS}$: free dyadic space $Ca^{2+}$ concentration.

Sarcoplasmic reticulum (SR) is a membrane structure within the muscle cell, whose main function is to store $Ca^{2+}$. Dyadic space is the region bounded by the T-tubule and SR. $Ca^{2+}$ ions are considered buffered if they are bound to negatively charged proteins (called buffers). Otherwise they are considered free. The action potential transmitted through the gap junction cause the current myocyte to depolarize. The depolarization opens the L-type (long-lasting) Ca channel located on the surface membrane. A small amount of $Ca^{2+}$ enters the myocyte due to potential, leading to a sharp increase of $Ca^{2+}$ in the dyadic space, which is a small region. This increase triggers the SR to release a large amount of $Ca^{2+}$ (calcium-induced calcium release) to enable the excitation-contraction. During diastole, calcium is removed from the cytoplasm through two ways. Ca is pumped (1) back into the SR and (2) out of the cell, mainly by the sodium-calcium exchange (NCX).

For the TTP model in svFSI, the following units have to be used: time in [ms], length [mm], amount of substance [mmol], voltage [mV]. Mass can be in [g].

Structures involved in $Ca^{2+}$ cycling[6].

Phenomenological Models

Phenomenological models are derived based on some observations of the full ion model [4]. Instead of following the transmembrane ionic currents, they use an oscillation system with a fast ($V$) and a slow ($r$) variable to mimic the behaviors of the action potentials. The oscillators, without considering the diffusion term, are modeled as

$$ \frac{\mathrm{d}V}{\mathrm{d}t}=f^{V}(V,r) $$ $$ \frac{\mathrm{d}r}{\mathrm{d}t}=f^{r}(V,r) $$

Note that this set of equations describes the electrophysiology in a single cardiac myocyte, and the choice of $f^V$ and $f^r$ will determine if this is a pacemaker cell or a non-pacemaker cell.

The FitzHugh-Nagumo (FN) model describes the pacemaker cells with

$$ f^V = c[V(V-\alpha)(1-V)-r] $$ $$ f^r = V-br+a $$

Action potential calculated from FitzHugh-Nagumo model.

The Aliev-Panfilov (AP) model describes the non-pacemaker cells with

$$ f^V = cV(V-\alpha)(1-V)-Vr $$

$$ f^r = \left( \gamma+\frac{\mu_1r}{\mu_2+V}\right)[-r-cV(V-b-1)] $$

Action potential calculated from Aliev-Panfilov model.

Note that $V$ and $t$ are non-dimensional values here. The following equations are used to recover the physiological action potential and time:

$$ \mathrm{FitzHugh-Nagumo model}: V^{fhn}=(65V-35)mV; ~~ t^{fhn} = (220t) ms $$

$$ \mathrm{Aliev-Panfilov model}: V^{ap}=(100V-80)mV; ~~ t^{ap} = (12.9t) ms $$

Another class of phenomenological models exists that include additional variables to account for intra-cellular calcium kinetics. One such model is the Bueno-Orovio (BO) model [5] that can serve as a trade-off between the complex TTP model and the simplified phenomenological AP model.

Activation Models Available in svFSI

The following table provides a summary of all the available electrophysiology models in svFSI.

Available cardiac electrophysiology models.
Electrophysiology Model Input Keyword
Aliev-Panfilov model[4] “ap”, “aliev-panfilov”
Fitzhugh-Nagumo model[4] “fn”, “fitzhugh-nagumo”
Bueno-Orovio-Cherry-Fenton model[5] “bo”, “bueno-orovio”
tenTusscher-Panfilov model[3] “ttp”, “tentusscher-panfilov”

Purkinje Network Plugin Tutorial

Introduction

The electrical activity in the heart tissue triggers muscle contraction and pumps blood into the systemic circulation. Under normal conditions, the electrical signal originates at the sinoatrial node located in the right atrium and reaches the atrioventricular node, which is the only electrical joint between the atria and the ventricles. [7] The bundle of His connects the atrioventricular node to a fast conducting network of fibers, called the Purkinje network, located beneath the inner-most layer of the heart wall. Purkinje cells are larger than the cardiomyocytes and conduct the excitation wave faster than any other cell on the heart tissue. [7]] The network not only synchronizes contraction between the left and the right ventricles but also allows the trajectory to follow a sequence beginning at the ventricular apex, spreading through the free-wall and eventually to the basal plane. Modeling the Purkinje network in cardiac electrophysiology simulations is therefore essential to achieve a realistic activation pattern and tissue contraction.

SimVascular provides a Purkinje plugin that can be used to generate Purkinje network on arbitrary patient-specific cardiac models. Our Purkinje plugin uses a fractal-tree-based method to generate the network [8] and provides a simple interface to adjust the parameters that control the network density and coverage. The Purkinje mesh is then exported in vtu format compatible with SimVascular/svFSI for the ensuing electrophysiology simulations.

In the following sections, we describe steps to download and install the Purkinje plugin, provide a workflow for creating Purkinje network using the plugin.

Plugin Installation

The installer for the Purkinje plugin can be found on the Simtk website. The installer installs the shared libraries, the setup.sh script and a Python script in,

/usr/local/sv/svplugins/SVDATE/Purkinje-Plugin/PLUGINDATE

Make sure that the SVDATE matches with the date of the installed SimVascular application. If there is a mismatch, you may rename the SVDATE folder to match with the installed SimVascular application. The setup.sh script sets the SVCUSTOMPLUGINS and SVPLUGINPATH environment variables. The SVCUSTOMPLUGINS environment variable defines the name of the plugin (e.g. orgsvguiqtpurkinjenetwork) and the SVPLUGINPATH environment variable defines the location of the plugin shared library (e.g. liborgsvguiqtpurkinjenetwork.so).

If successfully installed, the Purkinje Network tool icon automatically shows up on the main toolbar when the SimVascular application is launched.

Purkinje Plugin Workflow & Usage

The Purkinje plugin can be applied only on a SimVascular generated mesh and not on the model. Therefore, it is recommended that the user first create a mesh either by following the regular workflow (i.e., Image, Path, Segmentation, Model and Mesh) or by importing a solid model using the Models tool. In the example below, we use the latter approach and load the left ventricular endocardium directly into the SV project. We then mesh the volume before using the plugin to generate Purkinje network.

Create SV mesh

After creating an SV project, right-click on the Models tab to import a left ventricular solid model (.vtp) and run the Global Reinit command. During the import, extract faces on the model using the default settings.

Create a project or open an existing one.
Import .vtp surface file that you want to create a Purkinje network on.

Rename the surfaces in the imported solid model and select types from the drop down menu. Optionally assign colors to each surface to identify them easily. Identify the endocardium of the left ventricle and set it as a wall.

Rename the surfaces and set types to wall/cap where appropriate.

The imported left ventricular model may appear as shown below,

Imported left ventricular endocardium for creating the Purkinje network.

We now mesh the model before adding the Purkinje network tool. Create a new mesh under Meshes tab and open the SV Meshing window. Enter an appropriate max. edge size, deselect the Volume Meshing and run the mesher. Once the surface mesh is successfully generated, a window pops up with information about the mesh.

Create a mesh under Meshes tab.
Set mesh Global Max Edge Size and run the mesher.
Once the mesh is succesfully generated, the mesh stats window pops up.

Plugin workflow

Here we provide a brief overview of the workflow to generate Purkinje network. Once the mesh is created, activate the Purkinje plugin so that it displays the triangular surfaces of the mesh. Purkinje network can then be created using the following steps:

  1. Select a face on the mesh
  2. Select the network starting point on the face
  3. Set network generation parameters
  4. Generate the network
  5. Display the network
  6. Repeat the above steps until satisfied
  7. Export the network data

Below we provide a detailed description of the above steps with illustrations.

First, activate the Purkinje Network tool on the SimVascular toolbar that brings up the Purkinje Network tool panel. If the Purkinje Network window is already visible upon launching the SimVascular application, it is recommended to close the window and reactivate it only after an SV project has been created or loaded into the Data Manager.

Activate Purkinje tool by clicking the Purkinje Network icon shown on the toolbar.

A mesh face is selected by moving the mouse cursor over a face and pressing the s key. Select the endocardium face you want to generate the Purkinje network on. The selected face mesh is highlighed in yellow and its name is displayed in the Mesh Surface text box.

Select the surface to generate the Purkinje network on.

Provide a starting point of the Purkinje network by moving the mouse cursor to a node on the highlighted face mesh and pressing Ctrl(Cmd for OSX)+left-click. A red sphere marks the selected node and a red line segment showing the direction of the network is displayed. Adjust the starting point and direction by repeating the selection.

Select the starting point and direction of the Purkinje network.

Adjust the control parameters of the Purkinje network and generate by selecting the Create Network button. The network is displayed by checking the Show Network box. You can repeat the selection of the initial point and network creation until the results are satisfactory. A detailed description of the network control parameters is provided below.

Create the Purkinje network after adjusting parameters.
Purkinje network is overlaid on the mesh in the Display window.

The network is saved in .vtu format under the Purkinje-Network subfolder of the project directory. Additionally you can find the .txt files that contain the information of end nodes, position coordinates, and network connectivity under the same folder. Further details about the contents of each file is provided below.

The Purkinje network visualized in Paraview .

Output Description

The Purkinje Network tool writes a set of files in the Purkinje-Network subfolder within the SimVascular project directory. The files are prefixed with the name of the selected face. The files and their contents are described below:

  • FACENAME.vtp - triangular surface on which the network is generated.
  • FACENAME.vtu - network mesh represented as polylines of $n$ segments and $(n+1)$ nodes.
  • FACENAME_endnodes.txt - indices of nodes at the ends of network segments (i.e. not connected to other nodes).
  • FACENAME_ien.txt - network connectivity.
  • FACENAME_xyz.txt - network nodal coordinates.

Network Control Parameters

The Purkinje network is generated using a Python script. The inputs to the script include a triangular surface, the network starting point, a second point defining the direction of the first network branch, and the parameters used to control the shape, density and spread of the network [2]. These parameters are described below:

  • Starting point - initial node of the network.
  • Second point - defines the direction along which the network’s initial branch will grow.
  • Number of branch generations - denotes the maximum number of network branches generated from the initial node.
  • Average branch length - the length of each branch is calculated from a random normal distribution with a mean length, $L$ and a variance, $0.4 L^2$.
  • Branch angle - the angle between direction of the previous branch and a new branch.
  • Repulsive parameter - regulates the branch curvature (higher value leads to greater repulsion between the branches).
  • Branch segment length - an approximate length of the segments that compose one branch.

Known Issues

  1. SimVascular does not record that the plugin was added to a project. The plugin must be added to the project every time a project is opened. The project **Purkinje-Network** directory is saved between project sessions.
  2. The tool may not work if the Purkinje Network tool is added to a project before a mesh is created or loaded.
  3. The Purkinje Network tool does not know when the mesh changes. If the mesh is changed then the project must be saved and then reopened.
  4. There is no error reporting from the Python script. Users must check the console window for errors.
  5. The vtu mesh generated by the plugin cannot be read by svFSI as it is in ASCII format. A temporary solution is to export a regular appended/binary vtu from Paraview.

Examples

In this section, we will provide two examples of cardiac electrophysiology modeling. The first example is a widely used benchmark test, activation of a block of cardiac tissue [1]. The second example is the propagation of action potential in the Purkinje fibers. More examples are available on GitHub.

Activation of a Myocardial Block

The computational domain is a rectangular block, with dimensions of $3\times 7\times 20 mm$. A small cluster of cells located at one corner of the block (marked as S) will recive an initial stimulus, and the TTP model is used to model the activation of the rest of the domain. The entire simulation set-up along with results can be found on GitHub.

Computational domain of example 1.[reference]

Here is the breakdown of the input file for this study.


    # File svFSI.inp
    #----------------------------------------------------------------
    # General simulation parameters

    Continue previous simulation: f
    Number of spatial dimensions: 3
    Number of time steps: 100000
    Time step size: 0.01
    Spectral radius of infinite time step: 0.50
    Searched file name to trigger stop: STOP_SIM

    Save results to VTK format: 1
    Name prefix of saved VTK files: result
    Increment in saving VTK files: 100
    Start saving after time step: 1

    Increment in saving restart files: 100
    Convert BIN to VTK format: 0

    Verbose: 1
    Warning: 1
    Debug: 0

    #----------------------------------------------------------------
    # Mesh data
    Add mesh: msh {
       Mesh file path):   mesh/h0.1/mesh-complete.mesh.vtu
       Add face: X0 {
          Face file path: mesh/h0.1/mesh-surfaces/X0.vtp
       }
       Add face: X1 {
          Face file path: mesh/h0.1/mesh-surfaces/X1.vtp
       }
       Add face: Y0 {
          Face file path: mesh/h0.1/mesh-surfaces/Y0.vtp
       }
       Add face: Y1 {
          Face file path: mesh/h0.1/mesh-surfaces/Y1.vtp
       }
       Add face: Z0 {
          Face file path: mesh/h0.1/mesh-surfaces/Z0.vtp
       }
       Add face: Z1 {
          Face file path: mesh/h0.1/mesh-surfaces/Z1.vtp
       }

       Fiber direction: (1, 0, 0)

       # Here we also need to supply a domain information to separate
       # pacemaker cells and non-pacemaker cells. domain_info.dat
       # is a plaintext file store the domain ID of each element.
       # id == 1: non-pacemaker cells
       # id == 2: pacemaker cells (recivies stimulus)
       Domain file path: mesh/h0.1/domain_info.dat
    }

    #----------------------------------------------------------------
    # Equations
    Add equation: CEP {
       Coupled: 1
       Min iterations: 1
       Max iterations: 5
       Tolerance: 1e-6

       # Here domain id == 1 are non-pacemaker cells
       Domain: 1 {
          Electrophysiology model: TTP
          Conductivity (iso): 0.012571
          Conductivity (ani): 0.082715
          ODE solver: RK
       }

       # Here domain id == 2 are pacemaker cells
       Domain: 2 {
          Electrophysiology model: TTP
          Conductivity (iso): 0.012571
          Conductivity (ani): 0.082715
          # Stimulus to initiate the depolarization process
          Stimulus: Istim {
             Amplitude: -35.714
             Start time: 0.0
             Duration: 2.0
             Cycle length: 10000.0
          }
          ODE solver: RK
       }

       Output: Spatial {
          Action_potential: t
       }

       LS type: GMRES
       {
          Max iterations:      100
          Tolerance:           1D-6
          Krylov space dimension: 50
       }
    }
Activation of the cuboid.

Activation of a Purkinje Network

Here we perform electrophysiology simulation on a human biventricular Purkinje network. The initial stimulus is provided at the atrio-ventricular junction, and the TTP model is used to model for the activation of the individual Purkinje cells. The entire simulation set-up along with results can be found on GitHub.

A biventricular Purkinje network

Here is the breakdown of the input file for this study.


    #----------------------------------------------------------------
    # General simulation parameters

    Continue previous simulation: 1
    Number of spatial dimensions: 3
    Number of time steps: 20000
    Time step size: 0.05
    Spectral radius of infinite time step: 0.50
    Searched file name to trigger stop: STOP_SIM

    Save results to VTK format: 1
    Name prefix of saved VTK files: result
    Increment in saving VTK files: 50
    Start saving after time step: 1

    Increment in saving restart files: 1000
    Convert BIN to VTK format: 0

    Verbose: 1
    Warning: 0
    Debug: 0

    #----------------------------------------------------------------
    # Mesh data
    Add mesh: pfib_LV {
       Set mesh as fibers: t
       Mesh file path: mesh/pfib_LV.vtu
       Domain file path: mesh/domain_pfib_LV.dat
       Mesh scale factor: 10.0
    }

    Add mesh: pfib_RV {
       Set mesh as fibers: t
       Mesh file path: mesh/pfib_RV.vtu
       Domain file path: mesh/domain_pfib_RV.dat
       Mesh scale factor: 10.0
    }

    #----------------------------------------------------------------
    # Equations
    Add equation: CEP {
       Coupled: 1
       Min iterations: 1
       Max iterations: 3
       Tolerance: 1e-6

       Domain: 1 {
          Electrophysiology model: TTP
          Myocardial zone: pfib
          Conductivity (iso): 1.1
          ODE solver: RK
       }

       Domain: 2 {
          Electrophysiology model: TTP
          Myocardial zone: pfib
          Conductivity (iso): 1.1
          ODE solver: RK
          Stimulus: Istim {
             Amplitude: -52.0
             Start time: 0.0
             Duration: 2.0
             Cycle length: 1000.0
          }
       }

       Output: Spatial {
          Action_potential: t
       }

       LS type: GMRES
       {
          Max iterations:      100
          Tolerance:           1D-6
          Krylov space dimension: 50
       }
    }
Activation of the Purkinje network simulated using svFSI.

References

[1] Niederer, S. A., Kerfoot, E., Benson, A. P., Bernabeu, M. O., Bernus, O., Bradley, C., Cherry, E. M., Clayton, R., Fenton, F. H., Garny, A., Heidenreich, E., Land, S., Maleckar, M., Pathmanathan, P., Plank, G., Rodriguez, J. F., Roy, I., Sachse, F. B., Seemann, G., … Smith, N. P. (2011). Verification of cardiac tissue electrophysiology simulators using an N-version benchmark. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 369(1954), 4331–4351.

[2] ten Tusscher, K. H. W. J., Noble, D., Noble, P. J., & Panfilov, A. V. (2004). A model for human ventricular tissue. American Journal of Physiology-Heart and Circulatory Physiology, 286(4), H1573–H1589.

[3] ten Tusscher, K. H. W. J., & Panfilov, A. V. (2006). Alternans and spiral breakup in a human ventricular tissue model. American Journal of Physiology-Heart and Circulatory Physiology, 291(3), H1088–H1100.

[4] Göktepe, S., & Kuhl, E. (2009). Computational modeling of cardiac electrophysiology: A novel finite element approach. International Journal for Numerical Methods in Engineering, 79(2), 156–178.

[5] Bueno-Orovio, A., Cherry, E. M., & Fenton, F. H. (2008). Minimal model for human ventricular action potentials in tissue. Journal of Theoretical Biology, 253(3), 544–560.

[6] Bers, D. M. (2008). Calcium Cycling and Signaling in Cardiac Myocytes. Annual Review of Physiology, 70(1), 23–49.

[7] Dubin, D., (1996). Rapid Interpretation of EKG’s. Cover Publishing Company, Tampa, Florida.

[8] Sahli Costabal, F., Hurtado, D. E., & Kuhl, E. (2016). Generating Purkinje networks in the human heart. Journal of Biomechanics, 49(12), 2455–2465. https://doi.org/10.1016/j.jbiomech.2015.12.025






Cardiac Mechanics Modeling

Human heart is an intrinsic multiphysics system. It’s normal function requires synchronization and coupling of cardiac electrophysiology, tissue (myocardium) mechanics, and hemodynamics. Among the three physics, cardiac mechanics plays a central role, as it serves as the bridge between electrophysiology and blood flow. In a physiological setting, the myocardial cells (myocytes) in the heart tissue initiate and propagate an electrical impulse (action potential) that travels through the heart and drives the contraction of the cardiac muscle, which, in turn, propels the blood to sustain whole body circulation.

svFSI is a new multi-physics finite element solver designed for computational modeling of the whole heart dynamics. As the next generation underlining solver for the SimVascular software, svFSI is capable of simulating myocardial electrical activity using state-of-the-art cellular activation models, employ nonlinear hyperelastic material models to simulate tissue mechanics and activation coupling strategies, and perform large deformation fluid-structure interaction (FSI) to simulate the dynamics between heart tissue and blood flow. svFSI is compatible with several mesh formats and can read a wide range of 2D and 3D element types. The methods employed are widely popular among the scientific community and continue to evolve through research and development. The solver is parallelized using MPI for inter-process communication and was demonstrated to scale efficiently on large supercomputing clusters. ParMETIS is used for partitioning the computational domains. Simulation results are outputted into VTK format that can be easily visualized in the free software Paraview. In this document, we focus on demonstrating svFSI’s capability of cardiac mechanics modeling.

Solve Cardiac Mechanics with svFSI

Solving cardiac mechanics is a challenging task as the heart muscle is a complex fibrous structure that is predominantly incompressible [1], and undergoes large deformation during a cardiac cycle. In the following sections, we will introduce the workflow of modeling cardiac mechanics in svFSI using the benchmark example from Land et al. [2], i.e., the passive inflation of an idealized left ventricle.

Obtain svFSI

The source code of svFSI is publicly released through GitHub, and many new features are introduced frequently. Users are encouraged to download the most recent version and test it on their problems. svFSI is built through a CMake system and a short build guide is provided here. In addition to the source code, pre-built binary files are also available through SimTK. Here users can download svFSI binary for Ubuntu or MacOS.

The svFSI executable (either built from source or downloaded from SimTK) is not packaged into each SimVascular release, so users have to manually add its path to SimVascular GUI. The location is under Window -> Preference -> svFSI Simulation.

Set up path to **svFSI** solver in SimVascular.

Model Construction & Mesh Generation

SimVascular GUI provides a complete pipeline that includes model construction, mesh generation, problem configuration, input file generation, and numerical simulation. Although the model construction and mesh generation have been discussed extensively in other parts of the SimVascular user guides, we will go over them briefly in this subsection for the sake of completeness.

First, we create a SV Project called LV_Inflation. In cardiovascular modeling, computational models are usually generated from medical images. However, as we are using an idealized ventricular geometry here, we can directly import the model into the SV Data Manager by right clicking Models and selecting Import Solid Model. The solid model can be a surface mesh stored in either stl or vtp format.

Import solid model for mesh generation.

We could name the new model LV. When asked if you would like to extract faces from the model, click yes and accept the default separation angle (50 degree). Oftentimes, the imported model is not centered in the field of view. One can right click on the model and click Global Reinit to recenter the geometry. Double click the model name to open the SV Modeling window and rename each face corresponding to its physiological interpretation (base - LV basal plane; endo - endocardium; epi - epicardium) and set each of these face types as wall. You may also change the color of each face to distinguish them in the Display window.

Extract faces from the imported model.

Next we will generate an unstructured mesh from this idealized LV model. In the SV Data Manager, right click Meshes and select Create Mesh. Accepting the default options will create a mesh object named LV. Double click the mesh object to open the SV Meshing configuration window. The most import parameter here is Global Max Edge Size, which controls the size of the elements. We will set this value to 2.0 and run the mesher. A relatively coarse mesh will be generated. Users can reduce the Global Max Edge Size if a finer mesh is desired.

Meshing configuration.

To convert the generated mesh to the svFSI-ready format, right click the mesh object, LV, and select Export Mesh Complete.

Export the complete mesh.

Open the destination folder to find the mesh files contained in a folder called &ltmesh object name&gt-mesh-complete (LV-mesh-complete in this case). In this folder, the vtu file (mesh-complete.mesh.vtu) contains the volume mesh of the ventricle. The surface (boundary) mesh files could be found in the subfolder mesh-surfaces in vtp format. This set of vtk files defines the computational domain of this problem and can be visualized in Paraview.

Layout of the mesh folder.

Generation of Fiber Distribution

The heart wall is a composite of layers (or sheets) of parallel myocytes, which are the predominant fiber types. These fiber and sheet directions enable defining a local orthonormal coordinate system inside the cardiac muscle. This local coordinate system is crucial for using a structurally-based constitutive relation for the cardiac muscle [9]. For 3D problems, the users are required to prescribe fiber and sheet directions in the computational domain for certain constitutive relations. Users are provided an option to define a constant fiber direction through the following input directives,


# Fiber direction
Fiber direction: (1.0, 0.0, 0.0)

# Sheet direction
Fiber direction: (0.0, 1.0, 0.0)

svFSI also supports specifying distributed fiber and sheet direction generated by rule-based methods [3] through the following input commands,


# Fiber direction
Fiber direction file path: fibersLong.vtu

# Sheet direction
Fiber direction file path: fibersSheet.vtu

In the latter approach, the vtu file should have the local coordinate system stored as a vector with the name “FIB_DIR” at the centroid of each element. An example is provided here.

Input File

svFSI requires a plain-text input file to specify simulation parameters. An overview of the syntax could be found here. The SimVascular GUI currently supports limited input configurations. To access more advanced functions of svFSI, users are recommended to create their own input file by modifying existing templates. Below is a template of the input file for modeling passive inflation of a left ventricle,


# File svFSI.inp
#----------------------------------------------------------------
# General simulation parameters

Continue previous simulation: f
Number of spatial dimensions: 3
Number of time steps: 100
Time step size: 0.01
Spectral radius of infinite time step: 0.50
Searched file name to trigger stop: STOP_SIM

Save results to VTK format: 1
Name prefix of saved VTK files: result
Increment in saving VTK files: 1
Start saving after time step: 1

Increment in saving restart files: 100
Convert BIN to VTK format: 0

Verbose: 1
Warning: 1
Debug: 0

#----------------------------------------------------------------
# Mesh data including volume (vtu) and surface (vtp) meshes,
# domains and fiber distributions

Add mesh: msh {
   Mesh file path:    <path to mesh-complete folder>/mesh-complete.mesh.vtu
   Add face: endo {
      Face file path: <path to mesh-complete folder>/mesh-surfaces/endo.vtp
   }
   Add face: epi {
      Face file path: <path to mesh-complete folder>/mesh-surfaces/epi.vtp
   }
   Add face: base {
      Face file path: <path to mesh-complete folder>/mesh-surfaces/base.vtp
   }
   Fiber direction: (1.0, 0.0, 0.0)
   Fiber direction: (0.0, 1.0, 0.0)
}

#----------------------------------------------------------------
# Equations solved
# Here we use mixed formulation, ustruct, with stabilization
# displacement-based formulation can be invoked through struct
Add equation: ustruct {
   # Define min and max number of iterations, and convergence
   # tolerance for the nonlinear solver (Newton method)
   Coupled: 1
   Min iterations: 4
   Max iterations: 10
   Tolerance: 1e-6
   Use Taylor-Hood type basis: f

   # Define constitutive model and its parameters
   Constitutive model: Guccione {
      C:   1.0e4
      bf:  1.0
      bt:  1.0
      bfs: 1.0
   }

   # Define a small density value for quasi-steady (static)
   # simulation
   Density: 1e-3

   # Define Poisson ratio
   Poisson ratio: 0.5

#==================================================================================
#  This block is for struct; remove if ustruct is used
   # Define elasticity modulus for the volumetric part of the
   # strain energy function
   Elasticity modulus: 1.0
   # Penalty method to enforce incompressibility
   Dilational penalty model: ST91
   # Use a penalty parameter, if different from bulk modulus
   Penalty parameter: 1.0E6
#==================================================================================

#==================================================================================
#  This block is for ustruct; remove if struct is used
   # Define elasticity modulus to calculate tauM/tauC
   Elasticity modulus: 2.0e5
   # Stabilization paramters
   Momentum stabilization coefficient: 1e-3
   Continuity stabilization coefficient: 1e-3
#==================================================================================

   # Define variables for output
   Output: Spatial {
      Displacement: t
      Velocity: t
      Jacobian: t
   }

   # Linear solver parameter
   # Jacobi preconditioner (FSILS) is the default.
   # For ustruct, trilinos-ilu preconditioner is recommended.
   LS type: GMRES
   {
      Preconditioner:      FSILS
      Tolerance:           1D-6
      Max iterations:      1000
      Krylov space dimension: 50
   }

   # Apply zero displacement BC at the base
   Add BC: base {
      Type: Dir
      Value: 0.0
      Impose on state variable integral: t
   }

   # Add pressure load at the endo surface
   # Here we define a ramp function to linearly
   # increase pressure load from zero to the
   # desired value to increase solver stability.
   # Follower pressure load is set as the direction
   # of the applied load follows deformation.
   Add BC: endo {
      Type: Neu
      Time dependence: Unsteady
      Temporal values file path: <path to load.dat file>
      Ramp function: t
      Follower pressure load: t
   }
}

The applied load on the endocardial surface is provided in a file, load.dat, defining the ramp function. In this file, the first line specifies that there are two data points and the value will change linearly. The second and third line specify the time and the value at that time. Note that for a ramp function, the code expects data at two time points only. Should the simulation go beyond the last time stamp (t=1.0 in the current example), a constant value equal to the last extrema (1.0e4) will be applied.


2    1
0.0  0.0
1.0  1.0e4
# content of load.dat

Run Simulation

Though we can run the simulation through the GUI. It currently does not support some advanced features. Hence, we will use the aforementioned input file to directly run the simulation in terminals. The MPI run can be initiated through the following command


  mpiexec -np &ltnumber of MPI processes&gt  &ltpath to svFSI executable&gt  &ltpath to the input file&gt

Please remember to change the path to mesh folder and path to the load.dat in the input file according to your local configuration. svFSI will create a results directory called n-procs, where n is the number of MPI processes for the simulation. This directory will have vtu files that contain all the requested output fields and a log file called histor.dat.

svFSI Nonlinear Mechanics Solver

This section briefly explains the theory and implementation of the nonlinear solid dynamics solver in svFSI. Two types of formulations are provided in svFSI for modeling solid mechanics: STRUCT and uSTRUCT. STRUCT uses a pure displacement-based formulation of the balance laws (mass and momentum conservation principles), while uSTRUCT uses a mixed (velocity-pressure) formulation of the governing equations. The latter approach uses either a stabilized equal-order discretization for velocity and pressure function spaces, or the inf-sup-conforming Taylor-Hood-type finite element discretization.

Kinematics

During a cardiac cycle, the heart undergoes large deformations which can no longer be described by linear elasticity. Since the domain is continuously deforming, we introduce the concepts of a reference configuration (denoted $\Omega{0}$) and a current configuration (denoted $\Omega{t}$). The reference configuration is fixed for all times, and often refers to the initial geometry when $t=0$. We will use the vector $\mathbf{X}$ to denote the physical coordinates of the geometry in the reference configuration, and define the following relation:

$$ \begin{equation} \mathbf{x}(\mathbf{X}, t)=\mathbf{X}+\mathbf{u}(\mathbf{X}, t) \end{equation} $$

where $\mathbf{x}$ describes the physical coordinates of the geometry in the current configuration and $\mathbf{u}$ is the time-varying displacement vector field on $\Omega{0}$ which acts as a mapping between the reference and current configurations such that $\mathbf{u}: \Omega{0} \Rightarrow \Omega_{t} .$ We also define the following relationships and tensors, which are standard for describing nonlinear structural mechanics:

$$ \begin{equation} \begin{split} \ddot{\mathbf{u}}=\frac{d^{2} \mathbf{u}}{d t^{2}}, \quad \mathbf{F}=\frac{\partial \mathbf{x}}{\partial \mathbf{X}}, \quad & J=\operatorname{det}(\mathbf{F}) \ \mathbf{E}=\frac{1}{2}(\mathbf{C}-\mathbf{I}), \quad & \mathbf{C}=\mathbf{F}^{T} \mathbf{F} \end{split} \end{equation} $$

where $\ddot{\mathbf{u}}$ refers to the structural acceleration and the time derivative operator, $d^{2} / d t^{2}$, is applied on the reference configuration. $\mathbf{F}$ denotes the deformation gradient tensor and $J$ is the the determinant of $\mathbf{F}$ tensor that denotes the Jacobian of the transformation. $\mathbf{C}$ is the Cauchy-Green deformation tensor while $\mathbf{E}$ is the Green-Lagrange strain tensor.

STRUCT: displacement-based nonlinear elastodynamics

The STRUCT equation solves equations for nonlinear structural dynamics using finite element formulation. We start with the function spaces and weak form. We require that the trial and weighting functions satisfy their respective properties on the current domain. The strong form of the momentum balance is

$$ \begin{equation} \begin{split} \rho\ddot{\mathbf{u}} = \nabla{x}\cdot{\sigma} + \rho\mathbf{fb}, \ \sigma\cdot \mathbf{n} = \mathbf{h},~\mathrm{on}~\left(\Gamma{t}\right){h}. \end{split} \end{equation} $$

where $\sigma$ is the stress term in the current configuration.

The resulting weak form of the structural dynamics equations is

$$ \begin{equation} \int{\Omega{t}} \mathbf{w} \cdot \rho \ddot{\mathbf{u}} d \Omega+\int{\Omega{0}} \nabla{X} \mathbf{w}:(\mathbf{F} \mathbf{S}) d \Omega-\int{\Omega{t}} \mathbf{w} \cdot \rho \mathbf{fb} d \Omega-\int{\left(\Gamma{t}\right){h}} \mathbf{w} \cdot \mathbf{h} d \Gamma{h}=0 \end{equation} $$

The acceleration term (i.e. $\int{\Omega{t}} \mathbf{w} \cdot \rho \ddot{\mathbf{u}} d \Omega$), body forcing term (i.e. $\int{\Omega{t}} \mathbf{w} \cdot \mathbf{fb} d \Omega$), and the natural boundary condition term (i.e. $\int{\left(\Gamma{t}\right){h}} \mathbf{w} \cdot \mathbf{h} d \Gamma_{h}$) are all evaluated in the current configuration. These terms are commonly referred to as external work done on the structure by body forces and surface tractions. The remaining stress term in the equation is often referred to as internal work done on the structure by internal stresses, which we will treat differently here. We rewrite this in the reference configuration in terms of the deformation gradient and the second Piola-Kirchhoff stress tensor, which is commonly denoted as $\mathbf{S}$.

svFSI uses a splitting approach where the strain energy and the resulting second Piola-Kirchhoff stress tensor, $\mathbf{S}$, are decomposed into deviatoric (or isochoric, $\mathbf{S}{iso}$) and dilational (or volumetric, $\mathbf{S}{vol}$) components. The specific form of $\mathbf{S}$ will depend on the chosen constitutive model (isochoric and volumetric). More information on these Material models can be found in the literature [1]. It is noted that the symbol $E$ denotes the elastic modulus of the structure and is not to be confused with $\mathbf{E}$, which denotes the Green-Lagrange strain tensor, and $\nu$ represents Poisson’s ratio. Some key material parameters can then be defined as, $$ \begin{equation} \lambda=\frac{E v}{(1+v)(1-2 v)}, \quad \mu=\frac{E}{2(1+v)}, \quad \kappa=\lambda + \frac{2}{3} \mu \end{equation} $$

where, $\lambda$ and $\mu$ are the Lam&eacute’s first and second parameters, respectively, and $\kappa$ is the bulk modulus. The second Piola-Kirchhoff stress tensors for a few standard constitutive models are given as,

$$ \begin{equation} \begin{split} & \mathbf{S}^{StVK} & = 2 \mu \mathbf{E} + \lambda \operatorname{tr}(\mathbf{E}) \mathbf{I}, & \quad \textrm{St. Venant-Kirchhoff} \ & \mathbf{S}^{mStVK} & = \kappa \operatorname{ln}(J) \mathbf{C}^{-1} + \mu(\mathbf{C} - \mathbf{I}), & \quad \textrm{Modified St. Venant-Kirchhoff} \ & \mathbf{S_{iso}}^{nHK} & = \mu J^{-2/3} \left(\mathbf{I} - \frac{1}{3} \operatorname{tr}(\mathbf{C}) \mathbf{C}^{-1} \right), & \quad \textrm{Neo-Hookean} \end{split} \end{equation} $$

where $\mathbf{I}$ is the identity matrix. For the Neo-Hookean and other hyperelastic constitutive models, the $\mathbf{S}$ tensor is computed as $\mathbf{S} = \mathbf{S{iso}} + \mathbf{S{vol}}$, where $\mathbf{S_{vol}} = p J \mathbf{C}^{-1}$, and $p$ is the hydrostatic pressure computed based on the chosen dilational strain-energy function. See the section on Material models and the corresponding references for the available dilational penalty models in svFSI .

uSTRUCT: mixed formulation nonlinear elastodynamics

svFSI allows solving nonlinear elastodynamics using a mixed formulation where the structure’s velocity and pressure are the unknown degrees of freedom [4]. Two variants are available within this feature: (a) first, a novel variational multiscale stabilized (VMS) formulation that allows equal-order interpolation of velocity and pressure bases using a unified framework [4]; (b) second, using the classical inf-sup stable Taylor-Hood type elements where the velocity basis is derived from a function space that is one order higher relative to the pressure basis. In the displacement-based formulation, a hyperelastic material model assumes the existence of a Helmholtz free energy potential. However, uSTRUCT postulates hyperelasticity using Gibbs free energy potential [4] and takes the following additive decoupled form as, $$ \begin{equation} G(\overline{\mathbf{C}},p,T) = G{iso}(\overline{\mathbf{C}},T) + G{vol}(p,T) \end{equation} $$

where $G{vol}(p,T)=\int \rho^{-1}dp$, $\overline{\mathbf{C}}=J^{-2/3}\mathbf{C}$, $p$ is the pressure and $T$ is the temperature. It is worth mentioning that the free energy above is the specific free energy, i.e. the energy per unit mass. The free energy per unit volume is $G^R=\rho0G$, where $\rho_0$ is the density of the reference configuration. The Helmholtz free energy per unit volume can be obtained by a Legendre transformation of $-G^R$ as,

$$ \begin{equation} H^R(\overline{\mathbf{C}},J,T)=\sup_p\left(-pJ+G^R(\overline{\mathbf{C}},p,T) \right). \end{equation} $$

and due to the additive splitting of the Gibbs free energy, we have,

$$ \begin{equation} H^R(\overline{\mathbf{C}},J,T)=G{iso}^R(\overline{\mathbf{C}},T) + \supp\left(-pJ+G_{vol}^R(p,T) \right). \end{equation} $$

It is noted that Gibbs free energy naturally introduces pressure into the stress term. The governing equations for the motion of a continuum body in the current configuration are,

$$ \begin{equation} \begin{split} & \frac{\mathrm{d} \mathbf{u}}{\mathrm{d} t} - \mathbf{v} = \mathbf{0} \

& \beta(p) \frac{\mathrm{d} p}{\mathrm{d} t} + \nabla_x \cdot \mathbf{v} = 0  \\

& \rho(p) \frac{\mathrm{d} \mathbf{v}}{\mathrm{d} t} - \nabla_x \cdot \mathbf{\sigma_{dev}} + \nabla_x p - \rho(p) \mathbf{f_b} = \mathbf{0}.

\end{split} \end{equation} $$

The above system of equations represent the kinematic relation between displacement and velocity, balance of mass and linear momentum. $\sigma_{dev}$ is the deviatoric Cauchy stress, while $\rho$ and $\beta$ are the density and isothermal compressibility coefficient, respectively, defined as functions of pressure. The expressions for these quantities are given as follows,

$$ \begin{equation} \begin{split} & \rho(p) = \left( \frac{\mathrm{d} G{vol}}{\mathrm{d} p} \right)^{-1} \quad &, \quad & \beta(p) = \frac{1}{\rho} \frac{\mathrm{d} \rho}{\mathrm{d} p} = -\frac{\partial^2 G{vol}}{\partial p^2} / \frac{\partial G_{vol}}{\partial p} \

& \mathbf{\sigma_{dev}} = J^{-1} \bar{\mathbf{F}} \left( \mathbb{P}:\bar{\mathbf{S}} \right) \bar{\mathbf{F}}^T \quad &, \quad & \bar{\mathbf{S}} = 2 \frac{\partial G_{iso}^R}{\partial \bar{\mathbf{C}}} = 2 \frac{\partial (\rho_0 G_{iso})}{\partial \bar{\mathbf{C}}},

\end{split} \end{equation} $$

where $\mathbb{P} = \mathbb{I} - \frac{1}{3}\mathbf{C} \otimes \mathbf{C}^{-1}$ is the projection tensor, $\bar{\mathbf{F}} = J^{-1/3} \mathbf{F}$ and $\bar{\mathbf{C}} = J^{-2/3} \mathbf{C}$.

This mixed finite element problem is stabilized using variational multiscale method to allow using equal-order interpolating functions for velocity and pressure unknowns, employ linear elements and handle material incompressibility without suffering from locking issues. Defining an appropriate mixed function space, the stabilized weak form can then be written in the current configuration as,

$$ \begin{equation} \begin{split} & \mathbf{B}k := & \int{\Omegax} \mathbf{w}\mathbf{u} \cdot \left( \frac{\mathrm{d} \mathbf{u}}{\mathrm{d} t} - \mathbf{v} \right) \mathrm{d} \Omega_x = \mathbf{0} \

& \mathbf{B}_p := & \int_{\Omega_x} w_p \left( \beta(p) \frac{\mathrm{d} p}{\mathrm{d} t} + \nabla_x \cdot \mathbf{v} \right) \mathrm{d} \Omega_\mathbf{x} \\

& & + \sum_e \int_{\Omega_x^e} \tau_M^e \nabla_x w_p \cdot \left( \rho(p)\frac{\mathrm{d} \mathbf{v}}{\mathrm{d} t} - \nabla_x \cdot \mathbf{\sigma_{dev}} + \nabla_x p - \rho(p)\mathbf{f_b} \right) \mathrm{d} \Omega_x^e = 0 \\

& \mathbf{B}_m := & \int_{\Omega_x} \left( \mathbf{w}_\mathbf{v} \cdot \rho(p) \frac{\mathrm{d} \mathbf{v}}{\mathrm{d} t} + \nabla_x \mathbf{w}_\mathbf{v} : \mathbf{\sigma_{dev}} - \nabla_x \cdot \mathbf{w}_\mathbf{v} p - \mathbf{w}_\mathbf{v} \cdot \rho(p)\mathbf{f_b} \right) \mathrm{d} \Omega_x \\

 & & -\int_{\Gamma_x^h} \mathbf{w}_\mathbf{v} \cdot \mathbf{h} \mathrm{d} \Gamma_x + \sum_e \int_{\Omega_x^e} \tau_C \left(\nabla_x \cdot \mathbf{w}_\mathbf{v} \right) \left( \beta(p) \frac{\mathrm{d} p}{\mathrm{d} t} + \nabla_x \cdot \mathbf{v} \right) \mathrm{d} \Omega_x^e = 0.

\end{split}

\end{equation} $$

The stabilization parameters are chosen as,

$$ \begin{equation} \mathbf{\tau}M = \tauM\mathbf{I}{nd}, \quad \tauM = cm \frac{\Delta x^e}{c\rho}, \quad \tauC = c_c c\Delta x^e \rho \end{equation} $$

where, $\Delta x^e$ is the diameter of the circumscribing sphere of the tetrahedral element, $cm$ and $cc$ are two non-dimensional parameters, and $c$ is the maximum wave speed in the solid body. For compressible materials, $c$ is the bulk wave speed. Assuming isotropic small-strain linear elastic material, the bulk wave speed can be approximated as, $c=\sqrt{ \left( \lambda+2\mu \right) / \rho0}$, where $\lambda$ and $\mu$ are the Lam&eacute’s parameters. For incompressible materials, $c = \sqrt{\frac{\mu}{\rho0}}$ is the shear wave speed. Further details about the formulation, finite element discretization and time integration could be found in Liu et al. [4].

Material Models

Below is the list of available material constitutive models in svFSI :

Volumetric constitutive models for the struct/ustruct equations
Volumetric Model Input Keyword
Quadratic model “quad”, “Quad”, “quadratic”, “Quadratic”
Simo-Taylor91 model[5] “ST91”, “Simo-Taylor91”
Miehe94 model[6] “M94”, “Miehe94”
Isochoric constitutive models for the struct/ustruct equations
Isochoric Model Input Keyword
Saint Venant-Kirchhoff$^\dagger$ “stVK”, “stVenantKirchhoff”
modified Saint Venant-Kirchhoff$^\dagger$ “m-stVK”, “modified-stVK”, “modified-stVenantKirchhoff”
Neo-Hookean model “nHK”, “nHK91”, “neoHookean”, “neoHookeanSimo91”
Mooney-Rivlin model “MR”, “Mooney-Rivlin”
Holzapfel-Gasser-Ogden model [7] “HGO”
Guccione model [8] “Guccione”, “Gucci”
Holzapfel-Ogden model [9] “HO”, “Holzapfel”

$^\dagger$ These models are not available for ustruct.

References

[1] Holzapfel, G. A. (2002). Nonlinear solid mechanics: a continuum approach for engineering science. Wiley.

[2] Land, S., Gurev, V., Arens, S., Augustin, C. M., Baron, L., Blake, R., Bradley, C., Castro, S., Crozier, A., Favino, M., Fastl, T. E., Fritz, T., Gao, H., Gizzi, A., Griffith, B. E., Hurtado, D. E., Krause, R., Luo, X., Nash, M. P., … Niederer, S. A. (2015). Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 471 (2184), 20150641. https://doi.org/10.1098/rspa.2015.0641

[3] Bayer, J. D., Blake, R. C., Plank, G., & Trayanova, N. A. (2012). A Novel Rule-Based Algorithm for Assigning Myocardial Fiber Orientation to Computational Heart Models. Annals of Biomedical Engineering, 40(10), 2243–2254. https://doi.org/10.1007/s10439-012-0593-5

[4] Liu, J., & Marsden, A. L. (2018). A unified continuum and variational multiscale formulation for fluids, solids, and fluid–structure interaction. Computer Methods in Applied Mechanics and Engineering, 337, 549–597. https://doi.org/10.1016/J.CMA.2018.03.045

[5] Simo, J. C., & Taylor, R. L. (1991). Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms. Computer Methods in Applied Mechanics and Engineering, 85(3), 273–310. https://doi.org/10.1016/0045-7825(91)90100-K

[6] Miehe, C. (1994). Aspects of the formulation and finite element implementation of large strain isotropic elasticity. International Journal for Numerical Methods in Engineering, 37(12), 1981–2004. https://doi.org/10.1002/nme.1620371202

[7] Gasser, T. C., Ogden, R. W., & Holzapfel, G. A. (2006). Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Journal of The Royal Society Interface, 3(6), 15–35. https://doi.org/10.1098/rsif.2005.0073

[8] Guccione, J. M., McCulloch, A. D., & Waldman, L. K. (1991). Passive material properties of intact ventricular myocardium determined from a cylindrical model. Journal of Biomechanical Engineering, 113(February), 42–55. https://doi.org/10.1115/1.2894084

[9] Holzapfel, G. A, & Ogden, R. W. (2009). Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philosophical Transactions of the Royal Society Series A, 367(1902), 3445–3475. https://doi.org/10.1098/rsta.2009.0091

[10] Nolan, D. R., Gower, A. L., Destrade, M., Ogden, R. W., & McGarry, J. P. (2014). A robust anisotropic hyperelastic formulation for the modelling of soft tissue. Journal of the Mechanical Behavior of Biomedical Materials, 39, 48–60. https://doi.org/10.1016/J.JMBBM.2014.06.016