This tutorial was generated from the file projects/CaDiffusion/test/TestCaDiffusionLiteratePaper.hpp at revision r27328. Note that the code is given in full at the bottom of the page.
Ca^2+^ Channel Re-localization to Plasma-Membrane Microdomains Strengthens Activation of Ca^2+^-Dependent Nuclear Gene Expression
Code to accompany the paper Samanta et al. 2015.
Code Walkthrough
The following wiki page provides a walk-through of the Chaste code that was used to perform the simulations in this paper.
First we include some header files:
#include <cxxtest/TestSuite.h>
#include "GmshMeshReader.hpp"
#include "UblasIncludes.hpp"
//#include "SimpleLinearEllipticSolver.hpp"
#include "SimpleLinearParabolicSolver.hpp"
#include "TetrahedralMesh.hpp"
#include "BoundaryConditionsContainer.hpp"
#include "ConstBoundaryCondition.hpp"
#include "OutputFileHandler.hpp"
#include "RandomNumberGenerator.hpp"
#include "PetscSetupAndFinalize.hpp"
Set up a diffusion equation with a source term
d[Ca]/dt = D_Ca Laplacian([Ca]) + Q
[Ca] in units of uM D_Ca = 300 (nm)^2^/us integral of Q per ion channel’s worth of elements over which is to be applied = 2.5133e4 uM / us
template <unsigned SPACE_DIM>
class DiffusionEquationWithSourceTerm : public AbstractLinearParabolicPde<SPACE_DIM>
{
private:
const std::vector<c_vector<double, SPACE_DIM> >& mrChannelLocations;
/** The elements that are source elements */
std::vector<std::vector<unsigned> > mElementsEachChannel;
std::vector<double> mConcentrationSourcePerUnitVolume; // in units of uM / us
public:
DiffusionEquationWithSourceTerm(AbstractTetrahedralMesh<SPACE_DIM,SPACE_DIM>* pMesh,
const std::vector<c_vector<double, SPACE_DIM> >& rChannelLocations)
: AbstractLinearParabolicPde<SPACE_DIM>(),
mrChannelLocations(rChannelLocations)
{
// We're going to work out the total volume of the elements that are going to have a source term
// in them here.
std::vector<double> source_volumes(mrChannelLocations.size(), 0.0);
mConcentrationSourcePerUnitVolume.resize(mrChannelLocations.size());
mElementsEachChannel.resize(mrChannelLocations.size());
// Make a list of elements that we are going to say are source elements.
for (typename TetrahedralMesh<SPACE_DIM,SPACE_DIM>::ElementIterator elt_iter = pMesh->GetElementIteratorBegin();
elt_iter != pMesh->GetElementIteratorEnd();
++elt_iter)
{
// Decide whether this should be classed as a source element or not.
c_vector<double, SPACE_DIM> location = elt_iter->CalculateCentroid();
for (unsigned channel=0; channel < mrChannelLocations.size(); channel++)
{
if (norm_2(location-mrChannelLocations[channel]) <= 3.0) // If centroid is within 3nm of channel say it is a source.
{
//std::cout << "Source Element Recorded\n";
// Make a note of this element
mElementsEachChannel[channel].push_back(elt_iter->GetIndex());
// And add its volume to the total volume of source.
c_matrix<double, 3, 3> jacob;
double det;
elt_iter->CalculateJacobian(jacob, det);
source_volumes[channel] += elt_iter->GetVolume(det);
//std::cout << "Total source volume[" << channel << "] = " << source_volumes[channel] << "\n";
}
}
}
// This number is magical and should not be changed unless ion current changed.
const double total_source_required = 2.5133e4; // Fiddly conversion from ionic current in single channel.
for (unsigned channel=0; channel<mrChannelLocations.size(); channel++)
{
mConcentrationSourcePerUnitVolume[channel] = total_source_required / source_volumes[channel];
}
}
double ComputeSourceTerm(const ChastePoint<SPACE_DIM>& rPoint,
double u,
Element<SPACE_DIM,SPACE_DIM>* pElement)
{
for (unsigned channel=0; channel<mElementsEachChannel.size(); channel++)
{
if (std::find(mElementsEachChannel[channel].begin(), mElementsEachChannel[channel].end(), pElement->GetIndex())
!= mElementsEachChannel[channel].end())
{
return mConcentrationSourcePerUnitVolume[channel];
}
}
return 0.0;
}
The Diffusion constant for calcium is 300 um^2^ / s This is equivalent to 300 (nm)^2^ / us
c_matrix<double, SPACE_DIM, SPACE_DIM> ComputeDiffusionTerm(const ChastePoint<SPACE_DIM>& rPoint,
Element<SPACE_DIM,SPACE_DIM>* pElement=NULL)
{
const double diffusion_constant = 300; // Units : 3000 (nm)^2/us
return diffusion_constant*identity_matrix<double>(SPACE_DIM);
}
double ComputeDuDtCoefficientFunction(const ChastePoint<SPACE_DIM>& rPoint)
{
return 1.0;
}
};
Test class and method to look at Calcium diffusion
class TestCaDiffusion : public CxxTest::TestSuite
{
public:
We will solve du/dt = div(grad u) + u, in 3d, with boundary conditions Ca=0 on the boundary, and initial conditions Ca=0. *
Our units throughout this are: distance : nanometers nm time : microseconds us
void TestSolvingDiffusionEquationForCalcium() throw(Exception)
{
TetrahedralMesh<3,3> mesh;
Either a 3D Disc shaped thing, centred at (0,0,0), radius 10, and height 1.5. N.B. We’re in slightly odd units of 10 nm or 1e-8 m (!)
Note a coarse version of the mesh is provided to test simulations, but the ones in the paper were run on the refined version included here:
bool disc = true;
if (disc)
{
GmshMeshReader<3,3> gmsh_reader("projects/CaDiffusion/test/meshes/CaDiffusion.msh"); // for accurate solution
//GmshMeshReader<3,3> gmsh_reader("projects/CaDiffusion/test/meshes/CaDiffusionCoarse.msh"); // for quick estimate
mesh.ConstructFromMeshReader(gmsh_reader);
}
Or a square slab of membrane we construct on the fly
Create a 20 by 20 by 1.5 mesh in 3D, this time using the
[ConstructRegularSlabMesh](https://github.com/Chaste/trac_archive/wiki/Construct-Regular-Slab-Mesh)
method on the mesh. The first parameter is the cartesian space-step and the other three parameters are the width, height and depth of the mesh.
else
{
mesh.ConstructRegularSlabMesh(0.25, 20.0, 20.0, 1.5);
mesh.Translate(-10.0,-10.0,0.0); // Centre it in x-y plane at origin
}
mesh.Scale(10,10,10); // To get into units of nm, radius 100nm, height 15nm.
Create some ion channel locations
double z_location = 15.0; //nm - on the top / outer membrane.
std::vector<c_vector<double, 3u> > channel_locations;
// Create ion channels in a pentagon shape.
{
// 6.3 nm is the closest the channel pores could ever get from structures
// 9.6 nm is a value that was just published (Pemi et al. PNAS (2015) Nanoscale patterning of STIM1 and Orai1 during store-operated Ca entry)
// 47.5 nm is the mean nearest neighbour
// 88.5 nm is the mean between any two points (unlikely to be this spread).
double inter_channel_spacing = 9.6;
double pentagon_circumradius = (1.0/10.0)*(sqrt(50.0 + 10.0*sqrt(5.0)))*inter_channel_spacing;
std::cout << "Circumradius = " << pentagon_circumradius << std::endl;
// co-ordinates of a pentagon's vertices
double c1 = pentagon_circumradius*(cos(2.0*M_PI/5.0));
double c2 = pentagon_circumradius*(cos(M_PI/5.0));
double s1 = pentagon_circumradius*(sin(2.0*M_PI/5.0));
double s2 = pentagon_circumradius*(sin(4.0*M_PI/5.0));
c_vector<double, 3u> location;
location[0] = pentagon_circumradius;
location[1] = 0.0;
location[2] = z_location;
channel_locations.push_back(location);
location[0] = c1;
location[1] = s1;
location[2] = z_location;
channel_locations.push_back(location);
location[0] = -c2;
location[1] = s2;
location[2] = z_location;
channel_locations.push_back(location);
location[0] = -c2;
location[1] = -s2;
location[2] = z_location;
channel_locations.push_back(location);
location[0] = c1;
location[1] = -s1;
location[2] = z_location;
channel_locations.push_back(location);
}
Create the PDE object (defined above)
DiffusionEquationWithSourceTerm<3u> pde(&mesh, channel_locations);
Create a new boundary conditions container and specify u=0.0 on the boundary.
BoundaryConditionsContainer<3u,3u,1u> bcc; // Templated over element dim, space dim, problem dim.
ConstBoundaryCondition<3u>* p_bc_for_Ca = new ConstBoundaryCondition<3u>(0.0);
for (TetrahedralMesh<3u,3u>::BoundaryNodeIterator node_iter = mesh.GetBoundaryNodeIteratorBegin();
node_iter != mesh.GetBoundaryNodeIteratorEnd();
++node_iter)
{
double x = (*node_iter)->rGetLocation()[0];
double y = (*node_iter)->rGetLocation()[1];
//double z = (*node_iter)->rGetLocation()[2];
if (disc)
{
// If we are on the rim of the disc
if (x*x + y*y >= 100*100 - 1e-6)
{
bcc.AddDirichletBoundaryCondition(*node_iter, p_bc_for_Ca, 0);
}
}
else
{
// If we are on the edge of a slab
if (fabs(x+100) < 1e-6 || fabs(x-100) < 1e-6 || fabs(y+100)<1e-6 || fabs(y-100)<1e-6)
{
bcc.AddDirichletBoundaryCondition(*node_iter, p_bc_for_Ca, 0);
}
}
}
SimpleLinearParabolicSolver<3,3> solver(&mesh, &pde, &bcc);
For parabolic problems, initial conditions are also needed. The solver will expect a PETSc vector, where the i-th entry is the initial solution at node i, to be passed in. To create this PETSc
Vec
, we will use a helper function in the
[PetscTools](https://chaste.cs.ox.ac.uk/public-docs/classPetscTools.html)
class to create a
Vec
of size num_nodes, with each entry set to 0.0. Then we set the initial condition on the solver.
Vec initial_condition = PetscTools::CreateAndSetVec(mesh.GetNumNodes(), 0.0);
solver.SetInitialCondition(initial_condition);
Next define the start time, end time, and timestep, and set them.
double t_start = 0; // micro seconds
double t_end = 1; // micro seconds
double dt = 0.001; // micro seconds
solver.SetTimes(t_start, t_end);
solver.SetTimeStep(dt);
solver.SetOutputDirectoryAndPrefix("CaDiffusion/time_dependent","results");
solver.SetOutputToVtk(true);
Vec result = solver.Solve();
// Write a copy of the mesh to examine in a different format.
TrianglesMeshWriter<3,3> writer("CaDiffusion/mesh", "disc", false);
writer.WriteFilesUsingMesh(mesh);
All PETSc vectors should be destroyed when they are no longer needed.
PetscTools::Destroy(initial_condition);
PetscTools::Destroy(result);
}
};
Code
The full code is given below
File name TestCaDiffusionLiteratePaper.hpp
#include <cxxtest/TestSuite.h>
#include "GmshMeshReader.hpp"
#include "UblasIncludes.hpp"
//#include "SimpleLinearEllipticSolver.hpp"
#include "SimpleLinearParabolicSolver.hpp"
#include "TetrahedralMesh.hpp"
#include "BoundaryConditionsContainer.hpp"
#include "ConstBoundaryCondition.hpp"
#include "OutputFileHandler.hpp"
#include "RandomNumberGenerator.hpp"
#include "PetscSetupAndFinalize.hpp"
template <unsigned SPACE_DIM>
class DiffusionEquationWithSourceTerm : public AbstractLinearParabolicPde<SPACE_DIM>
{
private:
const std::vector<c_vector<double, SPACE_DIM> >& mrChannelLocations;
/** The elements that are source elements */
std::vector<std::vector<unsigned> > mElementsEachChannel;
std::vector<double> mConcentrationSourcePerUnitVolume; // in units of uM / us
public:
DiffusionEquationWithSourceTerm(AbstractTetrahedralMesh<SPACE_DIM,SPACE_DIM>* pMesh,
const std::vector<c_vector<double, SPACE_DIM> >& rChannelLocations)
: AbstractLinearParabolicPde<SPACE_DIM>(),
mrChannelLocations(rChannelLocations)
{
// We're going to work out the total volume of the elements that are going to have a source term
// in them here.
std::vector<double> source_volumes(mrChannelLocations.size(), 0.0);
mConcentrationSourcePerUnitVolume.resize(mrChannelLocations.size());
mElementsEachChannel.resize(mrChannelLocations.size());
// Make a list of elements that we are going to say are source elements.
for (typename TetrahedralMesh<SPACE_DIM,SPACE_DIM>::ElementIterator elt_iter = pMesh->GetElementIteratorBegin();
elt_iter != pMesh->GetElementIteratorEnd();
++elt_iter)
{
// Decide whether this should be classed as a source element or not.
c_vector<double, SPACE_DIM> location = elt_iter->CalculateCentroid();
for (unsigned channel=0; channel < mrChannelLocations.size(); channel++)
{
if (norm_2(location-mrChannelLocations[channel]) <= 3.0) // If centroid is within 3nm of channel say it is a source.
{
//std::cout << "Source Element Recorded\n";
// Make a note of this element
mElementsEachChannel[channel].push_back(elt_iter->GetIndex());
// And add its volume to the total volume of source.
c_matrix<double, 3, 3> jacob;
double det;
elt_iter->CalculateJacobian(jacob, det);
source_volumes[channel] += elt_iter->GetVolume(det);
//std::cout << "Total source volume[" << channel << "] = " << source_volumes[channel] << "\n";
}
}
}
// This number is magical and should not be changed unless ion current changed.
const double total_source_required = 2.5133e4; // Fiddly conversion from ionic current in single channel.
for (unsigned channel=0; channel<mrChannelLocations.size(); channel++)
{
mConcentrationSourcePerUnitVolume[channel] = total_source_required / source_volumes[channel];
}
}
double ComputeSourceTerm(const ChastePoint<SPACE_DIM>& rPoint,
double u,
Element<SPACE_DIM,SPACE_DIM>* pElement)
{
for (unsigned channel=0; channel<mElementsEachChannel.size(); channel++)
{
if (std::find(mElementsEachChannel[channel].begin(), mElementsEachChannel[channel].end(), pElement->GetIndex())
!= mElementsEachChannel[channel].end())
{
return mConcentrationSourcePerUnitVolume[channel];
}
}
return 0.0;
}
c_matrix<double, SPACE_DIM, SPACE_DIM> ComputeDiffusionTerm(const ChastePoint<SPACE_DIM>& rPoint,
Element<SPACE_DIM,SPACE_DIM>* pElement=NULL)
{
const double diffusion_constant = 300; // Units : 3000 (nm)^2/us
return diffusion_constant*identity_matrix<double>(SPACE_DIM);
}
double ComputeDuDtCoefficientFunction(const ChastePoint<SPACE_DIM>& rPoint)
{
return 1.0;
}
};
class TestCaDiffusion : public CxxTest::TestSuite
{
public:
void TestSolvingDiffusionEquationForCalcium() throw(Exception)
{
TetrahedralMesh<3,3> mesh;
bool disc = true;
if (disc)
{
GmshMeshReader<3,3> gmsh_reader("projects/CaDiffusion/test/meshes/CaDiffusion.msh"); // for accurate solution
//GmshMeshReader<3,3> gmsh_reader("projects/CaDiffusion/test/meshes/CaDiffusionCoarse.msh"); // for quick estimate
mesh.ConstructFromMeshReader(gmsh_reader);
}
else
{
mesh.ConstructRegularSlabMesh(0.25, 20.0, 20.0, 1.5);
mesh.Translate(-10.0,-10.0,0.0); // Centre it in x-y plane at origin
}
mesh.Scale(10,10,10); // To get into units of nm, radius 100nm, height 15nm.
double z_location = 15.0; //nm - on the top / outer membrane.
std::vector<c_vector<double, 3u> > channel_locations;
// Create ion channels in a pentagon shape.
{
// 6.3 nm is the closest the channel pores could ever get from structures
// 9.6 nm is a value that was just published (Pemi et al. PNAS (2015) Nanoscale patterning of STIM1 and Orai1 during store-operated Ca entry)
// 47.5 nm is the mean nearest neighbour
// 88.5 nm is the mean between any two points (unlikely to be this spread).
double inter_channel_spacing = 9.6;
double pentagon_circumradius = (1.0/10.0)*(sqrt(50.0 + 10.0*sqrt(5.0)))*inter_channel_spacing;
std::cout << "Circumradius = " << pentagon_circumradius << std::endl;
// co-ordinates of a pentagon's vertices
double c1 = pentagon_circumradius*(cos(2.0*M_PI/5.0));
double c2 = pentagon_circumradius*(cos(M_PI/5.0));
double s1 = pentagon_circumradius*(sin(2.0*M_PI/5.0));
double s2 = pentagon_circumradius*(sin(4.0*M_PI/5.0));
c_vector<double, 3u> location;
location[0] = pentagon_circumradius;
location[1] = 0.0;
location[2] = z_location;
channel_locations.push_back(location);
location[0] = c1;
location[1] = s1;
location[2] = z_location;
channel_locations.push_back(location);
location[0] = -c2;
location[1] = s2;
location[2] = z_location;
channel_locations.push_back(location);
location[0] = -c2;
location[1] = -s2;
location[2] = z_location;
channel_locations.push_back(location);
location[0] = c1;
location[1] = -s1;
location[2] = z_location;
channel_locations.push_back(location);
}
DiffusionEquationWithSourceTerm<3u> pde(&mesh, channel_locations);
BoundaryConditionsContainer<3u,3u,1u> bcc; // Templated over element dim, space dim, problem dim.
ConstBoundaryCondition<3u>* p_bc_for_Ca = new ConstBoundaryCondition<3u>(0.0);
for (TetrahedralMesh<3u,3u>::BoundaryNodeIterator node_iter = mesh.GetBoundaryNodeIteratorBegin();
node_iter != mesh.GetBoundaryNodeIteratorEnd();
++node_iter)
{
double x = (*node_iter)->rGetLocation()[0];
double y = (*node_iter)->rGetLocation()[1];
//double z = (*node_iter)->rGetLocation()[2];
if (disc)
{
// If we are on the rim of the disc
if (x*x + y*y >= 100*100 - 1e-6)
{
bcc.AddDirichletBoundaryCondition(*node_iter, p_bc_for_Ca, 0);
}
}
else
{
// If we are on the edge of a slab
if (fabs(x+100) < 1e-6 || fabs(x-100) < 1e-6 || fabs(y+100)<1e-6 || fabs(y-100)<1e-6)
{
bcc.AddDirichletBoundaryCondition(*node_iter, p_bc_for_Ca, 0);
}
}
}
SimpleLinearParabolicSolver<3,3> solver(&mesh, &pde, &bcc);
Vec initial_condition = PetscTools::CreateAndSetVec(mesh.GetNumNodes(), 0.0);
solver.SetInitialCondition(initial_condition);
double t_start = 0; // micro seconds
double t_end = 1; // micro seconds
double dt = 0.001; // micro seconds
solver.SetTimes(t_start, t_end);
solver.SetTimeStep(dt);
solver.SetOutputDirectoryAndPrefix("CaDiffusion/time_dependent","results");
solver.SetOutputToVtk(true);
Vec result = solver.Solve();
// Write a copy of the mesh to examine in a different format.
TrianglesMeshWriter<3,3> writer("CaDiffusion/mesh", "disc", false);
writer.WriteFilesUsingMesh(mesh);
PetscTools::Destroy(initial_condition);
PetscTools::Destroy(result);
}
};