| Chaste

This tutorial was generated from the file projects/CardiacEpVerification/test/TestConductionVelocityCaseStudyLiteratePaper.hpp at revision r20535. Note that the code is given in full at the bottom of the page.

Determining the accuracy of the conduction velocity in 1D

This is the main code for the first calculation verification case study. The simulation just involves a simple extension of a standard monodomain simulation.

The following are all standard includes:

#include <cxxtest/TestSuite.h>
#include "MonodomainProblem.hpp"
#include "DistributedTetrahedralMesh.hpp"
#include "AbstractCardiacCellFactory.hpp"
#include "LuoRudy1991.hpp"
#include "PetscSetupAndFinalize.hpp"

Define a simple cell factory which creates Luo-Rudy cells, and stimulates the given region.

class SimpleCellFactory : public AbstractCardiacCellFactory<1>
{
private:
    boost::shared_ptr<SimpleStimulus> mpStimulus;
    double mStimWidth;

public:
    SimpleCellFactory(double stimWidth)
        : AbstractCardiacCellFactory<1>(),
          mpStimulus(new SimpleStimulus(-5e5, 0.5)),
          mStimWidth(stimWidth)
    {
    }

    AbstractCardiacCell* CreateCardiacCellForTissueNode(unsigned nodeIndex)
    {
        double x = this->GetMesh()->GetNode(nodeIndex)->rGetLocation()[0];
        if (x<=mStimWidth)
        {
            return new CellLuoRudy1991FromCellML(mpSolver, mpStimulus);
        }
        else
        {
            return new CellLuoRudy1991FromCellML(mpSolver, mpZeroStimulus);
        }
    }
};

This class inherits from MonodomainProblem but does some extra work at the end of every (printing) timestep. (Note: the printing timestep will be set to be the same as the pde timestep).

class MonodomainProblemWithCvComputer1d : public MonodomainProblem<1>
{
private:
    double mQuarterNodeActivationTime;          // to be computed: activation time of node at x=0.25
    double mThreeQuarterNodeActivationTime;     // to be computed: activation time of node at x=0.75
    ReplicatableVector* mpVoltageLastTimestep;  // We will need to save the voltages at the last timestep
    double mLastTime;

public:
    MonodomainProblemWithCvComputer1d(AbstractCardiacCellFactory<1>* pCellFactory)
        : MonodomainProblem<1>(pCellFactory),
          mQuarterNodeActivationTime(DBL_MAX),
          mThreeQuarterNodeActivationTime(DBL_MAX),
          mpVoltageLastTimestep(NULL),
          mLastTime(DBL_MAX)
    {
    }

    ~MonodomainProblemWithCvComputer1d()
    {
        delete mpVoltageLastTimestep;
    }

At the end of every timestep, determine if the voltage for the nodes at x=0.25 and x=0.75 have just become positive. If so, use the current value of the voltage and the last value of the voltage, and linear interpolation over time, to get the time the voltage became positive.

    void OnEndOfTimestep(double time)
    {
        unsigned quarter_index = (this->mpMesh->GetNumNodes()-1)/4;
        unsigned three_quarter_index = 3*(this->mpMesh->GetNumNodes()-1)/4;

        double width = this->mpMesh->GetWidth(0/*ie x-direction*/);

        // a quick check that the nodes are in the expected places
        if((fabs(this->mpMesh->GetNode(quarter_index)->rGetLocation()[0]-0.25*width)>1e-12) || (fabs(this->mpMesh->GetNode(three_quarter_index)->rGetLocation()[0]-0.75*width)>1e-12))
        {
            std::cout << this->mpMesh->GetNode(quarter_index)->rGetLocation()[0] << " " << 0.25*width << " "
                      << this->mpMesh->GetNode(three_quarter_index)->rGetLocation()[0] << "\n";
            NEVER_REACHED;
        }

        ReplicatableVector voltage_repl(this->mSolution);
        if(mpVoltageLastTimestep==NULL)
        {
            mpVoltageLastTimestep = new ReplicatableVector(this->mSolution);
        }

        // determine activation times
        if(mQuarterNodeActivationTime==DBL_MAX && voltage_repl[quarter_index]>0)
        {
            double V2 = voltage_repl[quarter_index];
            double V1 = (*mpVoltageLastTimestep)[quarter_index];
            double t2 = time;
            double t1 = mLastTime;

            mQuarterNodeActivationTime = (V2*t1 - V1*t2)/(V2-V1);
        }

        if(mThreeQuarterNodeActivationTime==DBL_MAX && voltage_repl[three_quarter_index]>0)
        {
            double V2 = voltage_repl[three_quarter_index];
            double V1 = (*mpVoltageLastTimestep)[three_quarter_index];
            double t2 = time;
            double t1 = mLastTime;

            mThreeQuarterNodeActivationTime = (V2*t1 - V1*t2)/(V2-V1);
        }

        mLastTime = time;
    }

Get the conduction velocity from the activation times.

    double GetConductionVelocity()
    {
        unsigned quarter_index = (this->mpMesh->GetNumNodes()-1)/4;
        unsigned three_quarter_index = 3*(this->mpMesh->GetNumNodes()-1)/4;

        if( mQuarterNodeActivationTime==DBL_MAX || mThreeQuarterNodeActivationTime==DBL_MAX)
        {
            EXCEPTION("One of the nodes was not stimulated");
        }

        double dx = this->mpMesh->GetNode(three_quarter_index)->rGetLocation()[0] - this->mpMesh->GetNode(quarter_index)->rGetLocation()[0];
        return dx/(mThreeQuarterNodeActivationTime-mQuarterNodeActivationTime);
    }
};

class TestConductionVelocityCaseStudyLiteratePaper : public CxxTest::TestSuite
{
private:

The main simulation function:

    void Run(double parametersScaleFactor/*how much to scale h and dt*/, bool doTest=false /*see later*/)
    {

Define some initial parameters:

        double width = 1.0;      //cm
        double stim_width = 0.1; //cm
        double end_time = 10.0;  //ms

Define h and dt. Note h is proportional to dt.

        double init_h = 0.05; // cm, ie 500 um
        double h  = init_h*parametersScaleFactor;
        double dt = 0.01*parametersScaleFactor;

        double dt_ode = dt;
        double dt_pde = dt;
        double printing_dt = dt;

Run a standard monodomain simulation, except use our class MonodomainProblemWithCvComputer1d:

        DistributedTetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(h, width);

        HeartConfig::Instance()->SetSimulationDuration(end_time); //ms

        std::stringstream output_dir;
        output_dir << "CalculationVerification1d_" << parametersScaleFactor;
        HeartConfig::Instance()->SetOutputDirectory(output_dir.str());
        HeartConfig::Instance()->SetOutputFilenamePrefix("results");
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(dt_ode, dt_pde, printing_dt);

        SimpleCellFactory cell_factory(stim_width);

        MonodomainProblemWithCvComputer1d monodomain_problem( &cell_factory );
        monodomain_problem.SetMesh(&mesh);

        // only output at one node
        std::vector<unsigned> nodes_to_be_output;
        unsigned three_quarter_index = 3*(mesh.GetNumNodes()-1)/4;
        nodes_to_be_output.push_back(three_quarter_index);
        monodomain_problem.SetOutputNodes(nodes_to_be_output);

        //monodomain_problem.SetWriteInfo();
        monodomain_problem.Initialise();
        monodomain_problem.Solve();

Print results:

        std::cout << std::setprecision(9);
        double cv = monodomain_problem.GetConductionVelocity();
        std::cout << h << ", " << dt << ", " << cv << "\n";

If in ’testing mode’, which only applies if the coarsest mesh is being used, we do a quick test that nothing has changed:

        if(doTest)
        {
            if(parametersScaleFactor!=1.0)
            {
               EXCEPTION("Test with factor=1");
            }
            TS_ASSERT_DELTA(cv, 0.0736377, 1e-4);
        }
    }

The code which runs the above. The Richardson extrapolation of the results in done outside of Chaste. See the folder named ‘other’ in this project of a text file of the results.

public:
    void TestRunTest() throw (Exception)
    {
        Run(1.0, true);
    }

    void TestRun1dCv() throw (Exception)
    {
        unsigned num_sims = 5; // change to 9 for all paper results

        HeartConfig::Instance()->SetUseAbsoluteTolerance(1e-9);
        for(unsigned N=0; N<num_sims; N++)
        {
            double factor = 1.0/pow(2,N);
            Run(factor);
        }
    }
};

To obtain the results for the monodomain model problem QOI (Table 2), run TestMonodomain1d() in TestEpAgainstExactSolutionsLiteratePaper, after making a tiny edit to L2ErrorSquaredCalculator so that it only calculates the norm of the numerical solution, and ignores the exact solution part (see the comment in L2ErrorSquaredCalculator about Section 3.1).

Code

The full code is given below

File name `TestConductionVelocityCaseStudyLiteratePaper.hpp`

#include <cxxtest/TestSuite.h>
#include "MonodomainProblem.hpp"
#include "DistributedTetrahedralMesh.hpp"
#include "AbstractCardiacCellFactory.hpp"
#include "LuoRudy1991.hpp"
#include "PetscSetupAndFinalize.hpp"

class SimpleCellFactory : public AbstractCardiacCellFactory<1>
{
private:
    boost::shared_ptr<SimpleStimulus> mpStimulus;
    double mStimWidth;

public:
    SimpleCellFactory(double stimWidth)
        : AbstractCardiacCellFactory<1>(),
          mpStimulus(new SimpleStimulus(-5e5, 0.5)),
          mStimWidth(stimWidth)
    {
    }

    AbstractCardiacCell* CreateCardiacCellForTissueNode(unsigned nodeIndex)
    {
        double x = this->GetMesh()->GetNode(nodeIndex)->rGetLocation()[0];
        if (x<=mStimWidth)
        {
            return new CellLuoRudy1991FromCellML(mpSolver, mpStimulus);
        }
        else
        {
            return new CellLuoRudy1991FromCellML(mpSolver, mpZeroStimulus);
        }
    }
};

class MonodomainProblemWithCvComputer1d : public MonodomainProblem<1>
{
private:
    double mQuarterNodeActivationTime;          // to be computed: activation time of node at x=0.25
    double mThreeQuarterNodeActivationTime;     // to be computed: activation time of node at x=0.75
    ReplicatableVector* mpVoltageLastTimestep;  // We will need to save the voltages at the last timestep
    double mLastTime;

public:
    MonodomainProblemWithCvComputer1d(AbstractCardiacCellFactory<1>* pCellFactory)
        : MonodomainProblem<1>(pCellFactory),
          mQuarterNodeActivationTime(DBL_MAX),
          mThreeQuarterNodeActivationTime(DBL_MAX),
          mpVoltageLastTimestep(NULL),
          mLastTime(DBL_MAX)
    {
    }

    ~MonodomainProblemWithCvComputer1d()
    {
        delete mpVoltageLastTimestep;
    }

    void OnEndOfTimestep(double time)
    {
        unsigned quarter_index = (this->mpMesh->GetNumNodes()-1)/4;
        unsigned three_quarter_index = 3*(this->mpMesh->GetNumNodes()-1)/4;

        double width = this->mpMesh->GetWidth(0/*ie x-direction*/);

        // a quick check that the nodes are in the expected places
        if((fabs(this->mpMesh->GetNode(quarter_index)->rGetLocation()[0]-0.25*width)>1e-12) || (fabs(this->mpMesh->GetNode(three_quarter_index)->rGetLocation()[0]-0.75*width)>1e-12))
        {
            std::cout << this->mpMesh->GetNode(quarter_index)->rGetLocation()[0] << " " << 0.25*width << " "
                      << this->mpMesh->GetNode(three_quarter_index)->rGetLocation()[0] << "\n";
            NEVER_REACHED;
        }

        ReplicatableVector voltage_repl(this->mSolution);
        if(mpVoltageLastTimestep==NULL)
        {
            mpVoltageLastTimestep = new ReplicatableVector(this->mSolution);
        }

        // determine activation times
        if(mQuarterNodeActivationTime==DBL_MAX && voltage_repl[quarter_index]>0)
        {
            double V2 = voltage_repl[quarter_index];
            double V1 = (*mpVoltageLastTimestep)[quarter_index];
            double t2 = time;
            double t1 = mLastTime;

            mQuarterNodeActivationTime = (V2*t1 - V1*t2)/(V2-V1);
        }

        if(mThreeQuarterNodeActivationTime==DBL_MAX && voltage_repl[three_quarter_index]>0)
        {
            double V2 = voltage_repl[three_quarter_index];
            double V1 = (*mpVoltageLastTimestep)[three_quarter_index];
            double t2 = time;
            double t1 = mLastTime;

            mThreeQuarterNodeActivationTime = (V2*t1 - V1*t2)/(V2-V1);
        }

        mLastTime = time;
    }

    double GetConductionVelocity()
    {
        unsigned quarter_index = (this->mpMesh->GetNumNodes()-1)/4;
        unsigned three_quarter_index = 3*(this->mpMesh->GetNumNodes()-1)/4;

        if( mQuarterNodeActivationTime==DBL_MAX || mThreeQuarterNodeActivationTime==DBL_MAX)
        {
            EXCEPTION("One of the nodes was not stimulated");
        }

        double dx = this->mpMesh->GetNode(three_quarter_index)->rGetLocation()[0] - this->mpMesh->GetNode(quarter_index)->rGetLocation()[0];
        return dx/(mThreeQuarterNodeActivationTime-mQuarterNodeActivationTime);
    }
};

class TestConductionVelocityCaseStudyLiteratePaper : public CxxTest::TestSuite
{
private:
    void Run(double parametersScaleFactor/*how much to scale h and dt*/, bool doTest=false /*see later*/)
    {
        double width = 1.0;      //cm
        double stim_width = 0.1; //cm
        double end_time = 10.0;  //ms

        double init_h = 0.05; // cm, ie 500 um
        double h  = init_h*parametersScaleFactor;
        double dt = 0.01*parametersScaleFactor;

        double dt_ode = dt;
        double dt_pde = dt;
        double printing_dt = dt;

        DistributedTetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(h, width);

        HeartConfig::Instance()->SetSimulationDuration(end_time); //ms

        std::stringstream output_dir;
        output_dir << "CalculationVerification1d_" << parametersScaleFactor;
        HeartConfig::Instance()->SetOutputDirectory(output_dir.str());
        HeartConfig::Instance()->SetOutputFilenamePrefix("results");
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(dt_ode, dt_pde, printing_dt);

        SimpleCellFactory cell_factory(stim_width);

        MonodomainProblemWithCvComputer1d monodomain_problem( &cell_factory );
        monodomain_problem.SetMesh(&mesh);

        // only output at one node
        std::vector<unsigned> nodes_to_be_output;
        unsigned three_quarter_index = 3*(mesh.GetNumNodes()-1)/4;
        nodes_to_be_output.push_back(three_quarter_index);
        monodomain_problem.SetOutputNodes(nodes_to_be_output);

        //monodomain_problem.SetWriteInfo();
        monodomain_problem.Initialise();
        monodomain_problem.Solve();

        std::cout << std::setprecision(9);
        double cv = monodomain_problem.GetConductionVelocity();
        std::cout << h << ", " << dt << ", " << cv << "\n";

        if(doTest)
        {
            if(parametersScaleFactor!=1.0)
            {
               EXCEPTION("Test with factor=1");
            }
            TS_ASSERT_DELTA(cv, 0.0736377, 1e-4);
        }
    }

public:
    void TestRunTest() throw (Exception)
    {
        Run(1.0, true);
    }

    void TestRun1dCv() throw (Exception)
    {
        unsigned num_sims = 5; // change to 9 for all paper results

        HeartConfig::Instance()->SetUseAbsoluteTolerance(1e-9);
        for(unsigned N=0; N<num_sims; N++)
        {
            double factor = 1.0/pow(2,N);
            Run(factor);
        }
    }
};

Determining the accuracy of the conduction velocity in 1D#

Code#

File name TestConductionVelocityCaseStudyLiteratePaper.hpp#

Determining the accuracy of the conduction velocity in 1D

Code

File name `TestConductionVelocityCaseStudyLiteratePaper.hpp`