UserTutorials/RunningCryptSimulations

This tutorial is automatically generated from the file trunk/crypt/test/tutorial/TestRunningCryptSimulationsTutorial.hpp at revision r12206. Note that the code is given in full at the bottom of the page.

Examples showing how to run crypt simulations on periodic meshes with different cell-cycle models

Introduction

In this tutorial we show how Chaste can be used to simulate a cylindrical model of an intestinal crypt. Full details of the computational model can be found in the paper by van Leeuwen et al. (2009) [doi:10.1111/j.1365-2184.2009.00627.x].

The first thing to do is include the following header file, which allows us to use certain methods in our test. This header file must be included in any Chaste test.

#include <cxxtest/TestSuite.h>

Any test in which the GetIdentifier() method is used, even via the main cell_based code (through calls to AbstractCellPopulation output methods), must also include CheckpointArchiveTypes.hpp or CellBasedSimulationArchiver.hpp as the first Chaste header file.

#include "CheckpointArchiveTypes.hpp"

The next header file defines a helper class for generating cells for crypt simulations.

#include "CryptCellsGenerator.hpp"

The next two header files define two different types of cell-cycle model. In a FixedDurationGenerationBasedCellCycleModel, the duration of each phase of the cell cycle is fixed. In a WntCellCycleModel, the duration of a cell's G1 phase is determined by a system of nonlinear ODEs describing a cell's response to the local concentration of Wnt, a secreted cell–cell signalling molecule that is known to play a key role in cell proliferation in the crypt. In our crypt simulations, we impose a fixed gradient of Wnt up the axis of the crypt.

#include "FixedDurationGenerationBasedCellCycleModel.hpp"
#include "WntCellCycleModel.hpp"

The next header file defines a helper class for generating a suitable triangular mesh for the crypt simulation, such that the cell corresponding to each node is initially in mechanical equilibrium with its neighours and periodic boundary conditions are applied at the left- and right-hand sides of the mesh (hence the "cylindrical").

#include "CylindricalHoneycombMeshGenerator.hpp"

The next header file defines a CellPopulation class that uses a triangular mesh, and allows for the inclusion of 'ghost nodes'. These are nodes in the mesh that do not correspond to cells; instead they help ensure that a sensible Delaunay triangulation is generated at each timestep. This is because the triangulation algorithm requires a convex hull.

#include "MeshBasedCellPopulationWithGhostNodes.hpp"

The next header file defines a force law, based on a linear spring, for describing the mechanical interactions between neighbouring cells in the crypt.

#include "GeneralisedLinearSpringForce.hpp"

The next header file defines the class that simulates the evolution of a CellPopulation, specialized to deal with the cylindrical crypt geometry.

#include "CryptSimulation2d.hpp"

The next header file defines a Wnt singleton class, which (if used) deals with the imposed Wnt gradient in our crypt model. This affects cell proliferation in the case where we construct each cell with a WntCellCycleModel.

#include "WntConcentration.hpp"

The final header file defines a cell killer class, which implements sloughing of cells into the lumen once they reach the top of the crypt.

#include "SloughingCellKiller.hpp"

Next, we define the test class, which inherits from CxxTest::TestSuite and defines some test methods.

class TestRunningCryptSimulationsTutorial : public CxxTest::TestSuite
{
public:

Test 1: a basic crypt simulation

In the first test, we demonstrate how to create a crypt simulation using a cylindrical mesh, with each cell progressing through a fixed cell-cycle model, and sloughing enforced at the top of the crypt.

    void TestCryptWithFixedCellCycle() throw(Exception)
    {

As in all cell-based simulations, we must first set the start time.

        SimulationTime::Instance()->SetStartTime(0.0);

Next, we generate a mesh. The basic Chaste mesh is a TetrahedralMesh. To enforce periodicity at the left- and right-hand sides of the mesh, we use a subclass called Cylindrical2dMesh, which has extra methods for maintaining periodicity. To create a Cylindrical2dMesh, we can use a helper class called CylindricalHoneycombMeshGenerator. This generates a periodic honeycomb-shaped mesh, in which all nodes are equidistant to their neighbours. Here the first and second arguments define the size of the mesh - we have chosen a mesh that is 6 nodes (i.e. cells) wide, and 9 nodes high. The third argument indicates that we require a double layer of ghost nodes around the mesh (technically, just above and below the mesh, since it is periodic). We call GetCylindricalMesh() on the CylindricalHoneycombMeshGenerator to return our Cylindrical2dMesh, and call GetCellLocationIndices() to return a std::vector of indices of nodes in the mesh that correspond to real cells (as opposed to ghost nodes).

        CylindricalHoneycombMeshGenerator generator(6, 9, 2);
        Cylindrical2dMesh* p_mesh = generator.GetCylindricalMesh();
        std::vector<unsigned> location_indices = generator.GetCellLocationIndices();

Having created a mesh, we now create a std::vector of CellPtrs. To do this, we the CryptCellsGenerator helper class, which is templated over the type of cell-cycle model required (here FixedDurationGenerationBasedCellCycleModel) and the dimension. We create an empty vector of cells and pass this into the method Generate() along with the mesh. The third argument 'true' indicates that the cells should be assigned random birth times, to avoid synchronous division. The cells vector is populated once the method Generate() is called. Note that we only ever deal with shared pointers to cells, named CellPtrs.

        std::vector<CellPtr> cells;
        CryptCellsGenerator<FixedDurationGenerationBasedCellCycleModel> cells_generator;
        cells_generator.Generate(cells, p_mesh, location_indices, true);

Now we have a mesh, a set of cells to go with it, and a vector of node indices corresponding to real cells, we can create a CellPopulation object. In general, this class associates a collection of cells with a set of nodes or a mesh. For this test, because we have a mesh and ghost nodes, we use a particular type of cell population called a MeshBasedCellPopulationWithGhostNodes.

        MeshBasedCellPopulationWithGhostNodes<2> cell_population(*p_mesh, cells, location_indices);

Next we use the CellPopulation object to construct a CryptSimulation2d object, which will be used to simulate the crypt model.

        CryptSimulation2d simulator(cell_population);

We must set the output directory on the simulator (relative to "/tmp/<USER_NAME>/testoutput") and the end time (in hours).

        simulator.SetOutputDirectory("CryptTutorialFixedCellCycle");
        simulator.SetEndTime(1);

For longer simulations, we may not want to output the results every time step. In this case we can use the following method, to print results every 10 time steps instead. As the time step used by the simulator, is 30 seconds, this method will cause the simulator to print results every 5 minutes.

        simulator.SetSamplingTimestepMultiple(10);

Before running the simulation, we must add one or more force laws, which determine the mechanical behaviour of the cell population. For this test, we use a GeneralisedLinearSpringForce, which assumes that every cell experiences a force from each of its neighbours that can be represented as a linear overdamped spring.

        GeneralisedLinearSpringForce<2> linear_force;
        simulator.AddForce(&linear_force);

We also add a cell killer to the simulator. This object dictates under what conditions cells die. For this test, we use a SloughingCellKiller, which kills cells above a certain height (passed as an argument to the constructor).

        double crypt_height = 8.0;
        SloughingCellKiller<2> killer(&cell_population, crypt_height);
        simulator.AddCellKiller(&killer);

To run the simulation, we call Solve().

        simulator.Solve();

SimulationTime::Destroy() must be called at the end of the test. If not, when SimulationTime::Instance()->SetStartTime(0.0); is called at the beginning of the next test in this file, an assertion will be triggered.

        SimulationTime::Destroy();
    }

Finally, to visualize the results, we open a new terminal, cd to the Chaste directory, then cd to anim. Then we do: java Visualize2dCentreCells /tmp/$USER/testoutput/CryptTutorialFixedCellCycle/results_from_time_0. It may be necessary to do: javac Visualize2dCentreCells.java beforehand to create the java executable. Further details on visualization can be found on the Chaste wiki page For further details on visualization, see ChasteGuides/RunningCellBasedVisualization.

Test 2: a Wnt-dependent crypt simulation

The next test is very similar to Test 1, except that instead of using a fixed cell-cycle model, we use a Wnt-dependent cell-cycle model, with the Wnt concentration varying within the crypt in a predefined manner.

    void TestCryptWithWntCellCycle() throw(Exception)
    {

We first re-initialize time to zero and reseed the random number generator.

        SimulationTime::Instance()->SetStartTime(0.0);
        RandomNumberGenerator::Instance()->Reseed(0);

We then create a cylindrical mesh, and get the cell location indices, exactly as before.

        CylindricalHoneycombMeshGenerator generator(6, 9, 2);
        Cylindrical2dMesh* p_mesh = generator.GetCylindricalMesh();
        std::vector<unsigned> location_indices = generator.GetCellLocationIndices();

We create the cells, using the same method as before. Here, though, we use a WntCellCycleModel.

        std::vector<CellPtr> cells;
        CryptCellsGenerator<WntCellCycleModel> cells_generator;
        cells_generator.Generate(cells, p_mesh, location_indices, true);

We create the cell population, as before.

        MeshBasedCellPopulationWithGhostNodes<2> cell_population(*p_mesh, cells, location_indices);

We set the height of the crypt. As well as passing this variable into the sloughingCellKiller, we will pass it to the WntConcentration object (see below).

        double crypt_height = 8.0;

When using a WntCellCycleModel, we need a way of telling each cell what the Wnt concentration is at its location. To do this, we set up a WntConcentration object. Like SimulationTime, WntConcentration is a singleton class, so when instantiated it is accessible from anywhere in the code (and in particular, all cells and cell-cycle models can access it). We need to say what the profile of the Wnt concentation should be up the crypt: here, we say it is LINEAR (linear decreasing from 1 to 0 from the bottom of the crypt to the top). We also need to inform the WntConcentration of the cell population and the height of the crypt.

        WntConcentration<2>::Instance()->SetType(LINEAR);
        WntConcentration<2>::Instance()->SetCellPopulation(cell_population);
        WntConcentration<2>::Instance()->SetCryptLength(crypt_height);

Create a simulator as before (except setting a different output directory).

        CryptSimulation2d simulator(cell_population);
        simulator.SetOutputDirectory("CryptTutorialWntCellCycle");
        simulator.SetEndTime(1);

As before, we create a force law and cell killer and pass these objects to the simulator, then call Solve().

        GeneralisedLinearSpringForce<2> linear_force;
        simulator.AddForce(&linear_force);
        SloughingCellKiller<2> killer(&cell_population, crypt_height);
        simulator.AddCellKiller(&killer);
        simulator.Solve();

Finally, we must tidy up by destroying the SimulationTime and the WntConcentration singleton objects. This avoids memory leaks occurring.

        WntConcentration<2>::Destroy();
        SimulationTime::Destroy();
    }
};

The results of this test can be visualized as in Test 1, with the correct output directory.

Code

The full code is given below

File name `TestRunningCryptSimulationsTutorial.hpp`