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.

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();
    }
};

Cardiac Electro Mechanics

Cardiac electrophysiology: spiral wave

Cardiac electro-mechanics: a cardiac tissue wedge with varying fibre directions

Code overview#

Code#

File name TestElectroMechanicsLiteratePaper.hpp#

Code overview

Code

File name `TestElectroMechanicsLiteratePaper.hpp`