Cardiac electro-mechanics: a cardiac tissue wedge with varying fibre directions
This tutorial was generated from the file TestElectroMechanicsLiteratePaper.hpp. Note that the code is given in full at the bottom of the page.
![Electromechanics](/fig/paper-tutorials/composite.png)
On this wiki page we describe in detail the code that is used to run this example from the paper. This example is based on UserTutorials/CardiacElectroMechanics
The example shows how an electrical wave is followed by contraction along cardiac fibres. The fibres have different orientations in different parts of the tissue, resulting in a “twisting” motion.
Remember to run with build=GccOptNative
for speed.
e.g.
scons build=GccOptNative test_suite=projects/project_Plos2013/test/TestElectroMechanicsLiteratePaper.hpp
This example uses only files from the core Chaste code.
There is a special cmgui visualization script that loads both mechanics and voltage solution onto the
same mesh, which is in the same folder as this file: LoadElectroMechanicsSimulation.com
Code overview
We first include some header files which define the classes we will use.
#include <cxxtest/TestSuite.h>
#include "PlaneStimulusCellFactory.hpp"
#include "LuoRudy1991.hpp"
#include "NonlinearElasticityTools.hpp"
#include "CardiacElectroMechanicsProblem.hpp"
#include "FileFinder.hpp"
#include "VtkMeshWriter.hpp"
#include "PetscSetupAndFinalize.hpp"
class TestElectroMechanicsLiteratePaper : public CxxTest::TestSuite
{
public:
The following code is the test
itself, we use the scons
/ cxx-test
framework to run simulations, as it
provides a handy way to do all the necessary linking and library building.
void TestTwistingCube() throw(Exception)
{
Set the length of the simulation, long enough for some contraction to occur
HeartConfig::Instance()->SetSimulationDuration(36.0); //ms
This class sets up cells with the correct stimuli to stimulate just those cells on the x=0 surface at t=0.
We define the strength of stimulus that we need.
PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 3u> cell_factory(-1000*1000);
For electro-mechanics we need to set up two meshes of 1mm by 1mm by 1mm We have a fine electrics mesh and a coarser mechanics one.
TetrahedralMesh<3u,3u> electrics_mesh;
electrics_mesh.ConstructRegularSlabMesh(0.01, 0.1, 0.1, 0.1);
QuadraticMesh<3> mechanics_mesh(0.02, 0.1, 0.1, 0.1);
We are going to fix some nodes on Z (dimension ‘2’ indexed from 0) so we first need to identify these
std::vector<unsigned> fixed_nodes
= NonlinearElasticityTools<3u>::GetNodesByComponentValue(mechanics_mesh, 2, 0.0);
We now define the electro-mechanics problem. This includes specifying the contraction model and the material properties, as well as those nodes that are fixed and the mechanics timestep to use.
Mechanics happens over a longer timescale than electrophysiology, so we can use larger space and time steps for this aspect.
ElectroMechanicsProblemDefinition<3u> problem_defn(mechanics_mesh);
problem_defn.SetContractionModel(KERCHOFFS2003,1.0);
problem_defn.SetUseDefaultCardiacMaterialLaw(COMPRESSIBLE);
problem_defn.SetZeroDisplacementNodes(fixed_nodes);
problem_defn.SetMechanicsSolveTimestep(1.0); // ms
std::string output_directory = "Plos2013_ElectroMechanics";
std::string fibre_file_name = "5by5by5_fibres.ortho";
This is how to generate a fibre file for this mesh (this could be generated once, stored and re-loaded, as indeed we do below) We use a Streeter-style formula here.
Usually the fibres for a scientific problem would be determined by e.g. DTMRI and stored with the mesh files.
{
OutputFileHandler handler(output_directory + "Fibres");
out_stream p_file = handler.OpenOutputFile(fibre_file_name);
*p_file << mechanics_mesh.GetNumElements() << "\n"; // first line is number of entries
std::vector<c_vector<double,3u> > fibre_directions;
for(unsigned i=0; i<mechanics_mesh.GetNumElements(); i++)
{
double big_x = mechanics_mesh.GetElement(i)->CalculateCentroid()(0);
double theta = M_PI/3.0 - 10*big_x*2*M_PI/3.0; // 60 degrees when X=0, -60 when X=0.1;
c_vector<double,3u> fibre_direction;
fibre_direction[0] = 0;
fibre_direction[1] = cos(theta);
fibre_direction[2] = sin(theta);
*p_file << fibre_direction[0] << " " << fibre_direction[1] << " " << fibre_direction[2] // first three entries are fibre direction
<< " 0 " << -fibre_direction[2] << " " << fibre_direction[1] // next three are sheet direction
<< " 1 0 0\n"; // then normal to sheet direction
fibre_directions.push_back(fibre_direction);
}
p_file->close();
We only compile the following if VTK is installed and set up. This is optional - it is only used here for visualizing the fibre directions as in the paper figure. The simulation will run without it.
#ifdef CHASTE_VTK
VtkMeshWriter<3u,3u> mesh_writer(output_directory+ "Fibres", "mesh", false);
mesh_writer.AddCellData("Fibre Directions", fibre_directions);
mesh_writer.WriteFilesUsingMesh(mechanics_mesh);
#endif // CHASTE_VTK
}
Load up the file we just wrote to use as fibre directions for this mechanics problem.
FileFinder fibre_file_finder(output_directory + "Fibres/" + fibre_file_name, RelativeTo::ChasteTestOutput);
problem_defn.SetVariableFibreSheetDirectionsFile(fibre_file_finder, false);
Set up and solve the full cardiac electro-mechanics problem.
CardiacElectroMechanicsProblem<3u> problem(COMPRESSIBLE,
MONODOMAIN,
&electrics_mesh,
&mechanics_mesh,
&cell_factory,
&problem_defn,
output_directory);
Run the simulation
problem.Solve();
Report where time was spent to std::cout.
MechanicsEventHandler::Headings();
MechanicsEventHandler::Report();
}
};
Code
The full code is given below
File name TestElectroMechanicsLiteratePaper.hpp
#include <cxxtest/TestSuite.h>
#include "PlaneStimulusCellFactory.hpp"
#include "LuoRudy1991.hpp"
#include "NonlinearElasticityTools.hpp"
#include "CardiacElectroMechanicsProblem.hpp"
#include "FileFinder.hpp"
#include "VtkMeshWriter.hpp"
#include "PetscSetupAndFinalize.hpp"
class TestElectroMechanicsLiteratePaper : public CxxTest::TestSuite
{
public:
void TestTwistingCube() throw(Exception)
{
HeartConfig::Instance()->SetSimulationDuration(36.0); //ms
PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 3u> cell_factory(-1000*1000);
TetrahedralMesh<3u,3u> electrics_mesh;
electrics_mesh.ConstructRegularSlabMesh(0.01, 0.1, 0.1, 0.1);
QuadraticMesh<3> mechanics_mesh(0.02, 0.1, 0.1, 0.1);
std::vector<unsigned> fixed_nodes
= NonlinearElasticityTools<3u>::GetNodesByComponentValue(mechanics_mesh, 2, 0.0);
ElectroMechanicsProblemDefinition<3u> problem_defn(mechanics_mesh);
problem_defn.SetContractionModel(KERCHOFFS2003,1.0);
problem_defn.SetUseDefaultCardiacMaterialLaw(COMPRESSIBLE);
problem_defn.SetZeroDisplacementNodes(fixed_nodes);
problem_defn.SetMechanicsSolveTimestep(1.0); // ms
std::string output_directory = "Plos2013_ElectroMechanics";
std::string fibre_file_name = "5by5by5_fibres.ortho";
{
OutputFileHandler handler(output_directory + "Fibres");
out_stream p_file = handler.OpenOutputFile(fibre_file_name);
*p_file << mechanics_mesh.GetNumElements() << "\n"; // first line is number of entries
std::vector<c_vector<double,3u> > fibre_directions;
for(unsigned i=0; i<mechanics_mesh.GetNumElements(); i++)
{
double big_x = mechanics_mesh.GetElement(i)->CalculateCentroid()(0);
double theta = M_PI/3.0 - 10*big_x*2*M_PI/3.0; // 60 degrees when X=0, -60 when X=0.1;
c_vector<double,3u> fibre_direction;
fibre_direction[0] = 0;
fibre_direction[1] = cos(theta);
fibre_direction[2] = sin(theta);
*p_file << fibre_direction[0] << " " << fibre_direction[1] << " " << fibre_direction[2] // first three entries are fibre direction
<< " 0 " << -fibre_direction[2] << " " << fibre_direction[1] // next three are sheet direction
<< " 1 0 0\n"; // then normal to sheet direction
fibre_directions.push_back(fibre_direction);
}
p_file->close();
#ifdef CHASTE_VTK
VtkMeshWriter<3u,3u> mesh_writer(output_directory+ "Fibres", "mesh", false);
mesh_writer.AddCellData("Fibre Directions", fibre_directions);
mesh_writer.WriteFilesUsingMesh(mechanics_mesh);
#endif // CHASTE_VTK
}
FileFinder fibre_file_finder(output_directory + "Fibres/" + fibre_file_name, RelativeTo::ChasteTestOutput);
problem_defn.SetVariableFibreSheetDirectionsFile(fibre_file_finder, false);
CardiacElectroMechanicsProblem<3u> problem(COMPRESSIBLE,
MONODOMAIN,
&electrics_mesh,
&mechanics_mesh,
&cell_factory,
&problem_defn,
output_directory);
problem.Solve();
MechanicsEventHandler::Headings();
MechanicsEventHandler::Report();
}
};