PDESolvers in CompuCell3D¶
One of the most important and time consuming parts of the CC3D simulation is to solve all sorts of Partial Differential Equations which describe behavior of certain simulation objects (usually chemical fields). Most of the CC3D PDE solvers solve PDE with diffusive terms. Because we are dealing with moving boundary condition problems it was easiest and probably most practical to use explicit scheme of Finite Difference method. Most of CC3D PDE solvers run on multi core architectures and we also have GPU solvers which run and high performance GPU’s and they also provide biggest speedups in terms of performance. Because CC3D solvers were implemented at different CC3D life cycle and often in response to particular user requests, CC3DML specification may differ from solver to solver. However, the basic structure of CC3DML PDE solver code follows similar pattern .
- FlexibleDiffusionSolver
- Instabilities of the Forward Euler Method
- DiffusionSolverFE
- AdvectionDiffusionSolver.
- FastDiffusionSolver2D
- KernelDiffusionSolver
- ReactionDiffusionSolver
- Steady State diffusion solver
- Fluctuation Compensator Solver Add-On
FlexibleDiffusionSolver¶
This steppable is one of the basic and most important modules in CompuCell3D simulations.
Tip
Starting from version 3.6.2
we developed DiffusionSolverFE
which eliminates several inconveniences of FlexibleDiffusionSolver
.
As the name suggests it is responsible for solving diffusion equation but in addition to this it also handles chemical secretion which maybe thought of as being part of general diffusion equation:
where k
is a decay constant of concentration c
and D
is the
diffusion constant. The term called secretion
has the meaning as
described below.
Example syntax for FlexibleDiffusionSolverFE looks as follows:
<Steppable Type="FlexibleDiffusionSolverFE">
<AutoscaleDiffusion/>
<DiffusionField Name="FGF8">
<DiffusionData>
<FieldName>FGF8</FieldName>
<DiffusionConstant>0.1</DiffusionConstant>
<DecayConstant>0.002</DecayConstant>
<ExtraTimesPerMCS>5</ExtraTimesPerMCS>
<DeltaT>0.1</DeltaT>
<DeltaX>1.0</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
<InitialConcentrationExpression>x*y
</InitialConcentrationExpression>
</DiffusionData>
<SecretionData>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
<BoundaryConditions>
<Plane Axis="X">
<ConstantValue PlanePosition="Min" Value="10.0"/>
<ConstantValue PlanePosition="Max" Value="10.0"/>
</Plane>
<Plane Axis="Y">
<ConstantDerivative PlanePosition="Min" Value="10.0"/>
<ConstantDerivative PlanePosition="Max" Value="10.0"/>
</Plane>
</BoundaryConditions>
</DiffusionField>
<DiffusionField Name="FGF">
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>0.02</DiffusionConstant>
<DecayConstant>0.001</DecayConstant>
<DeltaT>0.01</DeltaT>
<DeltaX>0.1</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<SecretionOnContact Type="Medium"
SecreteOnContactWith="Amoeba">0.1</SecretionOnContact>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
</DiffusionField>
</Steppable>
We define sections that describe a field on which the steppable is to
operate. In our case we declare just two diffusion fields - FGF8
and FGF
.
Warning
When you want to solve more than one diffusive field using given diffusion
you should declare those multiple fiuelds within a single definition of the diffusion solver. You cannot
include two diffusion solver definitions - one with e.g. FGF8
and another with FGF
.
In other words you can only define e.g. FlexibleDiffusionSolverFE
once per simulation.
You can however use FlexibleDiffusionSolverFE for e.g. FGF
and DiffusionSolverFE
for e.g. FGF8
Inside the diffusion field we specify sections describing diffusion and
secretion. Let’s take a look at DiffusionData
section first:
<DiffusionField Name="FGF8">
<DiffusionData>
<FieldName>FGF8</FieldName>
<DiffusionConstant>0.1</DiffusionConstant>
<DecayConstant>0.002</DecayConstant>
<ExtraTimesPerMCS>5</ExtraTimesPerMCS>
<DeltaT>0.1</DeltaT>
<DeltaX>1.0</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
<InitialConcentrationExpression>x*y
</InitialConcentrationExpression>
</DiffusionData>
We give a name (FGF8
) to the diffusion field – this is required as we
will refer to this field in other modules.
Next we specify diffusion constant and decay constant.
Warning
We use Forward Euler Method to solve these equations. This is not a stable method for solving diffusion equation and we do not perform stability checks. If you enter too high diffusion constant for example you may end up with unstable (wrong) solution. Always test your parameters to make sure you are not in the unstable region.
We may also specify cells that will not participate in the diffusion.
You do it using <DoNotDiffuseTo>
tag. In this example we do not let any FGF
diffuse
into Bacteria
cells. You may of course use as many as necessary
<DoNotDiffuseTo>
tags. To prevent decay of a chemical in certain cells
we use syntax:
<DoNotDecayIn>Medium</DoNotDecayIn>
In addition to diffusion parameters we may specify how secretion should
proceed. SecretionData
section contains all the necessary information to
tell CompuCell how to handle secretion. Let’s study the example:
<SecretionData>
<SecretionOnContact Type="Medium" SecreteOnContactWith="Amoeba">0.1</SecretionOnContact>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
Here we have a definition two major secretion modes. Line:
<Secretion Type="Amoeba">0.1</Secretion>
ensures that every cell of type Amoeba
will get 0.1
increase in
concentration every MCS. Line:
<SecretionOnContact Type="Medium" SecreteOnContactWith="Amoeba">0.1</SecretionOnContact>
means that cells of type Medium
will get additional 0.1
increase in
concentration but only when they touch cell of type Amoeba
. This mode of
secretion is called SecretionOnContact
.
We can also see new CC3DML tags <DeltaT>
and <DeltaX>
. Their values
determine the correspondence between MCS and actual time and between
lattice spacing and actual spacing size. In this example for the first
diffusion field one MCS corresponds to 0.``1 units of actual time and
lattice spacing is equal ``1
unit of actual length. What is happening here
is that the diffusion constant gets multiplied by:
DeltaT/(DeltaX * DeltaX)
provided the decay constant is set to 0. If the decay constant is not zero DeltaT appears additionally in the term (in the explicit numerical approximation of the diffusion equation solution) containing decay constant so in this case it is more than simple diffusion constant rescaling.
DeltaT
and DeltaX
settings are closely related to ExtraTimesPerMCS
setting which allows calling of diffusion (and only diffusion) more than
once per MCS. The number of extra calls per MCS is specified by the user
on a per-field basis using ExtraTimesPerMCS
tag.
Warning
When using ExtraTimesPerMCS
secretion functions will
called only once per MCS. This is different than using PDESolverCaller
where entire module is called multiple times (this include diffusion and
secretion for all fields).
Tip
We recommend that you stay away from redefining DeltaX
and
DeltaT
and assume that your diffusion/decay coefficients are expressed
in units of pixel (distance) and MCS (time). This way when you assign
physical time and distance units to MCS and pixels you can easily
obtain diffusion and decay constants. DeltaX
and DeltaT
introduce
unnecessary complications.
The AutoscaleDiffusion
tag tells CC3D to automatically rescale diffusion
constant when switching between square and hex lattices. In previous
versions of CC3D such scaling had to be done manually to ensure that
solutions diffusion of equation on different lattices match. Here we
introduced for user convenience a simple tag that does rescaling
automatically. The rescaling factor comes from the fact that the
discretization of the divergence term in the diffusion equation has
factors such as unit lengths, using surface are and pixel/voxel volume
in it. On a square lattice all those values have numerical value of 1.0
.
On hex lattice, and for that matter of non-square latticeses, only
pixel/voxel volume has numerical value of 1
. All other quantities have
values different than 1.0
which causes the necessity to rescale
diffusion constant. The detail of the hex lattice derivation will be
presented in the “Introduction to Hexagonal Lattices in CompuCell3D” -
http://www.compucell3d.org/BinDoc/cc3d_binaries/Manuals/HexagonalLattice.pdf .
The FlexibleDiffusionSolver
is also capable of solving simple coupled
diffusion type PDE of the form:
where \(m_c\), \(m_f\), \(n_c\) , \(n_f\), \(p_c\), \(p_d\) are coupling coefficients. To code the above equations in CC3DML syntax you need to use the following syntax:
<Steppable Type="FlexibleDiffusionSolverFE">
<DiffusionField>
<DiffusionData>
<FieldName>c</FieldName>
<DiffusionConstant>0.1</DiffusionConstant>
<DecayConstant>0.002</DecayConstant>
<CouplingTerm InteractingFieldName=”d” CouplingCoefficent=”0.1”/>
<CouplingTerm InteractingFieldName=”f” CouplingCoefficent=”0.2”/>
<DeltaT>0.1</DeltaT>
<DeltaX>1.0</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
</DiffusionField>
<DiffusionField>
<DiffusionData>
<FieldName>d</FieldName>
<DiffusionConstant>0.02</DiffusionConstant>
<DecayConstant>0.001</DecayConstant>
<CouplingTerm InteractingFieldName=”c” CouplingCoefficent=”-0.1”/>
<CouplingTerm InteractingFieldName=”f” CouplingCoefficent=”-0.2”/>
<DeltaT>0.01</DeltaT>
<DeltaX>0.1</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
</DiffusionField>
<DiffusionField>
<DiffusionData>
<FieldName>f</FieldName>
<DiffusionConstant>0.02</DiffusionConstant>
<DecayConstant>0.001</DecayConstant>
<CouplingTerm InteractingFieldName=”c” CouplingCoefficent=”-0.2”/>
<CouplingTerm InteractingFieldName=”d” CouplingCoefficent=”0.2”/>
<DeltaT>0.01</DeltaT>
<DeltaX>0.1</DeltaX>
<DoNotDiffuseTo>Bacteria</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<Secretion Type="Amoeba">0.1</Secretion>
</SecretionData>
</DiffusionField>
</Steppable>
As one can see the only addition that is required to couple diffusion equations has simple syntax:
<CouplingTerm InteractingFieldName=”c” CouplingCoefficent=”-0.1”/>
<CouplingTerm InteractingFieldName=”f” CouplingCoefficent=”-0.2”/>
Instabilities of the Forward Euler Method¶
Most of the PDE solvers in CC3D use Forward Euler explicit numerical scheme. This method is unstable for large diffusion constant. As a matter of fact using D=0.25 with pulse initial condition will lead to instabilities in 2D. To deal with this you would normally use implicit solvers however due to moving boundary conditions that we have to deal with in CC3D simulations, memory requirements, performance and the fact that most diffusion constants encountered in biology are quite low (unfortunately this is not for all chemicals e.g. oxygen ) we decided to use explicit scheme. If you have to use large diffusion constants with explicit solvers you need to do rescaling:
- Set D, Δt, Δx according to your model
- If
you will need to to call solver multiple times per MCS.
- Set
<ExtraTimesPerMCS>
to N-1 where:
and
Initial Conditions¶
We can specify initial concentration as a function of x
, y
, z
coordinates using <InitialConcentrationExpression>
tag we use muParser
syntax to type the expression. Alternatively we may use
ConcentrationFileName
tag to specify a text file that contains values of
concentration for every pixel:
<ConcentrationFileName>initialConcentration2D.txt</ConcentrationFileName>
The value of concentration of the last pixel read for a given cell
becomes an overall value of concentration for a cell. That is if cell
has, say 8
pixels, and you specify different concentration at every
pixel, then cell concentration will be the last one read from the file.
Concentration file format is as follows:
*x y z c*
where x
, y
, z
, denote coordinate of the pixel. c
is the value of the
concentration.
Example:
0 0 0 1.2
0 0 1 1.4
...
The initial concentration can also be input from the Python script (typically in the start function of the steppable) but often it is more convenient to type one line of the CC3DML script than few lines in Python.
Boundary Conditions¶
All standard solvers (Flexible
, Fast
, and ReactionDiffusion
) by default
use the same boundary conditions as the GGH simulation (and those are
specified in the Potts section of the CC3DML script). Users can,
however, override those defaults and use customized boundary conditions
for each field individually. Currently CompuCell3D supports the
following boundary conditions for the diffusing fields: periodic,
constant value (Dirichlet) and constant derivative (von Neumann). To
specify custom boundary condition we include <BoundaryCondition> section
inside <DiffusionField>
tags.
The <BoundaryCondition>
section describes boundary conditions along
particular axes. For example:
<Plane Axis="X">
<ConstantValue PlanePosition="Min" Value="10.0"/>
<ConstantValue PlanePosition="Max" Value="10.0"/>
</Plane>
specifies boundary conditions along the x
axis. They are Dirichlet-type
boundary conditions. PlanePosition='Min"
denotes plane parallel to yz
plane passing through x=0
. Similarly PlanePosition="Min"
denotes plane
parallel to yz
plane passing through x=fieldDimX-1
where fieldDimX
is x
dimension of the lattice.
By analogy we specify constant derivative boundary conditions:
<Plane Axis="Y">
<ConstantDerivative PlanePosition="Min" Value="10.0"/>
<ConstantDerivative PlanePosition="Max" Value="10.0"/>
</Plane>
We can also mix types of boundary conditions along single axis:
<Plane Axis="Y">
<ConstantDerivative PlanePosition="Min" Value="10.0"/>
<ConstantValue PlanePosition="Max" Value="0.0"/>
</Plane>
Here in the xz
plane at y=0
we have von Neumann boundary conditions but
at y=fieldFimY-1
we have dirichlet boundary condition.
To specify periodic boundary conditions along, say, x
axis we use the
following syntax:
<Plane Axis="X">
<Periodic/>
</Plane>
Notice, that <Periodic>
boundary condition specification applies to both
“ends” of the axis i.e. we cannot have periodic boundary conditions at
x=0
and constant derivative at x=fieldDimX-1
.
DiffusionSolverFE¶
DiffusionSolverFE is new solver in 3.6.2 and is intended to fully
replace FlexibleDiffusionSolverFE
. It eliminates several limitations
and inconveniences of FlexibleDiffusionSolverFE
and provides new
features such as GPU implementation or cell type dependent
diffusion/decay coefficients. In addition it also eliminates the need
to rescale diffusion/decay/secretion constants. It checks stability
condition of the PDE and then rescales appropriately all coefficients
and computes how many extra times per MCS the solver has to be called. It
makes those extra calls automatically.
Warning
One of the key differences between FlexibleDiffusionSolverFE
and
DiffusionSolverFE
is the way in which secretion is treated. In
FlexibleDiffusionSolverFE
all secretion amount is done once followed by
possibly multiple diffusion calls to diffusion (to avoid numerical
instabilities). In DiffusionSolverFE
the default mode of operation is
such that multiple secretion and diffusion calls are interleaved.
This means that instead of secreting full amount for a given MCS and
diffusing it, the DiffusionSolverFE
secretes substance gradually so that
there is equal amount of secretion before each call of the diffusion.
One can change this behavior by adding <DoNotScaleSecretion/>
to
definition of the diffusion solver e.g.
With such definition the DiffusionSolverFE
will behave like
FlexibleDiffusionSolverFE
as far as computation.
Note
DiffusionSolverFE
autoscales diffusion discretization
depending on the lattice so that <AutoscaleDiffusion/>
we used in
FlexibleDiffusionSolverFE
is unnecessary.
This may result in slow performance so users have to be aware that those
extra calls to the solver may be the cause.
Typical syntax for the DiffusionSolverFE
may look like example below:
<Steppable Type="DiffusionSolverFE">
<DiffusionField Name="ATTR">
<DiffusionData>
<FieldName>ATTR</FieldName>
<GlobalDiffusionConstant>0.1</GlobalDiffusionConstant>
<GlobalDecayConstant>5e-05</GlobalDecayConstant>
<DiffusionCoefficient CellType="Red">0.0</DiffusionCoefficient>
</DiffusionData>
<SecretionData>
<Secretion Type="Bacterium">100</Secretion>
</SecretionData>
<BoundaryConditions>
<Plane Axis="X">
<Periodic/>
</Plane>
<Plane Axis="Y">
<Periodic/>
</Plane>
</BoundaryConditions>
</DiffusionField>
</Steppable>
The syntax resembles the syntax for FlexibleDiffusionSolverFE
. We
specify global diffusion constant by using <GlobalDiffusionConstant>
tag. This specifies diffusion coefficient which applies to entire region
of the simulation. We can override this specification for regions
occupied by certain cell types by using the following syntax:
<DiffusionCoefficient CellType="Red">0.0</DiffusionCoefficient>
Similar principles apply to decay constant and we use
<GlobalDecayConstant>
tag to specify global decay coefficient and
<DecayCoefficient CellType="Red">0.0</DecayCoefficient>
to override global definition for regions occupied by Red cells.
We do not support <DeltaX>
, <DeltaT>
or <ExtraTimesPerMCS>
tags.
Note
DiffusionSolverFE
autoscales diffusion discretization
depending on the lattice so that <AutoscaleDiffusion/>
we used in
FlexibleDiffusionSolverFE
is unnecessary.
Running DiffusionSolver on GPU¶
To run DiffusionSolverFE
on GPU all we have to do (besides having OpenCL
compatible GPU and correct drives installed) to replace first line of
solver specification:
<Steppable Type="DiffusionSolverFE">
with
<Steppable Type="DiffusionSolverFE_OpenCL">
Note
Depending on your computer hardware you may or may not be able to take advantage of GPU capabilities.
AdvectionDiffusionSolver.¶
Note
This is an experimental module and was not fully curated.
AdvectionDiffusionSolver
solves advection diffusion equation on a cell field as
opposed to grid. Of course, the inaccuracies are bigger than in the case
of PDE being solved on the grid but on the other hand solving the PDE on
a cell field means that we associate concentration with a given cell (not
just with a lattice point). This means that as cells move so does the
concentration. In other words we get advection for free. The
mathematical treatment of this kind of approximation was described in
Phys. Rev. E 72, 041909 (2005) paper by D.Dan et al.
The equation solved by this steppable is of the type:
where \(c\) denotes concentration , \(D\) is diffusion constant, \(k\) decay constant, \(\vec{\nu}\) is velocity field.
In addition to just solving advection-diffusion equation this module
allows users to specify secretion rates of the cells as well as
different secretion modes. The syntax for this module follows same
pattern as FlexibleDiffusionSolverFE
.
Example syntax:
<Steppable Type="AdvectionDiffusionSolverFE">
<DiffusionField Name="FGF">
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>0.05</DiffusionConstant>
<DecayConstant>0.003</DecayConstant>
<ConcentrationFileName>flowFieldConcentration2D.txt</ConcentrationFileName>
<DoNotDiffuseTo>Wall</DoNotDiffuseTo>
</DiffusionData>
<SecretionData>
<Secretion Type="Fluid">0.5</Secretion>
<SecretionOnContact Type="Fluid"
<SecreteOnContactWith="Wall">0.3</SecretionOnContact>
</SecretionData>
</DiffusionField>
</Steppable>
FastDiffusionSolver2D¶
Note
The current implementation of DiffusionSolverFE
may actually be faster than this legacy module. So
before commiting to it please verify the performance of teh DiffusionSolverFE
vs FastDiffusionSolverFE
FastDiffusionSolver2DFE
module is a simplified version of the
FlexibleDiffusionSolverFE
steppable. It runs several times faster that
flexible solver but lacks some of its features. Typical syntax is shown
below:
<Steppable Type="FastDiffusionSolver2DFE">
<DiffusionField Name="FGF">
<DiffusionData>
<UseBoxWatcher/>
<FieldName>FGF</FieldName>
<DiffusionConstant>0.010</DiffusionConstant>
<DecayConstant>0.003</DecayConstant>
<ExtraTimesPerMCS>2</ExtraTimesPerMCS>
<DoNotDecayIn>Wall</DoNotDecay>
<ConcentrationFileName>Demos/diffusion/diffusion_2D_fast_box.pulse.txt
</ConcentrationFileName>
</DiffusionData>
</DiffusionField>
</Steppable>
In particular, for fast solver you cannot specify cells into which diffusion is prohibited. However, you may specify cell types where diffusant decay is prohibited
For explanation on how <ExtraTimesPerMCS>
tag works works see section on
FlexibleDiffusionSolverFE
.
KernelDiffusionSolver¶
This diffusion solver has the advantage over previous solvers that it
can handle large diffusion constants. It is also stable. However, it
does not accept options like <DoNotDiffuseTo>
or <DoNotDecayIn>
. It also
requires periodic boundary conditions.
Simply put KernelDiffusionSolver solves diffusion equation:
with fixed, periodic boundary conditions on the edges of the lattice.
This is different from FlexibleDiffusionSolverFE
where the boundary
conditions evolve. You also need to choose a proper Kernel range (K)
according to the value of diffusion constant. Usually when \(K^2e^{-K^2/{4D}}\)
is small (this is the main part of the
integrand), the approximation converges to the exact value.
The syntax for this solver is as follows:
<Steppable Type="KernelDiffusionSolver">
<DiffusionField Name="FGF">
<Kernel>4</Kernel>
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>1.0</DiffusionConstant>
<DecayConstant>0.000</DecayConstant>
<ConcentrationFileName>
Demos/diffusion/diffusion_2D.pulse.txt
</ConcentrationFileName>
</DiffusionData>
</DiffusionField>
</Steppable>
Inside <DiffusionField>
tag one may also use option <CoarseGrainFactor>
. For example:
<Steppable Type="KernelDiffusionSolver">
<DiffusionField Name="FGF">
<Kernel>4</Kernel>
<CoarseGrainFactor>2</CoarseGrainFactor>
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>1.0</DiffusionConstant>
<DecayConstant>0.000</DecayConstant>
<ConcentrationFileName>
Demos/diffusion/diffusion_2D.pulse.txt
</ConcentrationFileName>
</DiffusionData>
</DiffusionField>
</Steppable>
ReactionDiffusionSolver¶
The reaction diffusion solver solves the following system of N reaction diffusion equations:
where W
denotes cell type
Let’s consider a simple example of such system:
It can be coded as follows:
<Steppable Type="ReactionDiffusionSolverFE">
<AutoscaleDiffusion/>
<DiffusionField Name="F">
<DiffusionData>
<FieldName>F</FieldName>
<DiffusionConstant>0.010</DiffusionConstant>
<ConcentrationFileName>
Demos/diffusion/diffusion_2D.pulse.txt
</ConcentrationFileName>
<AdditionalTerm>-0.01*H</AdditionalTerm>
</DiffusionData>
</DiffusionField>
<DiffusionField Name="H">
<DiffusionData>
<FieldName>H</FieldName>
<DiffusionConstant>0.0</DiffusionConstant>
<AdditionalTerm>0.01*F</AdditionalTerm>
</DiffusionData>
</DiffusionField>
</Steppable>
Notice how we implement functions f
from the general system of
reaction diffusion equations. We simply use <AdditionalTerm>
tag and
there we type arithmetic expression involving field names (tags
<FieldName>
). In addition to this we may include in those expression
word CellType
. For example:
<AdditionalTerm>0.01*F*CellType</AdditionalTerm>
This means that function f
will depend also on CellType
. CellType
holds the value of the type of the cell at particular location - x
, y
, z
- of the lattice. The inclusion of the cell type might be useful if you
want to use additional terms which may change depending of the cell
type. Then all you have to do is to either use if statements inside
<AdditionalTerm>
or form equivalent mathematical expression using
functions allowed by muParser
: http://muparser.sourceforge.net/mup_features.html#idDef2
For example, let’s assume that additional term for second equation is the following:
In such a case additional term would be coded as follows:
<AdditionalTerm>CellType==1 ? 0.01*F : 0.15*F</AdditionalTerm>
Notice that we have used here, so called ternary operator which might be familiar to you from other programing languages such as C or C++ and is equivalent to`` if-then-els``e statement
The syntax of the ternary (aka if-then-else
statement) is as follows:
condition ? expression if condition is true : expression if condition false
Warning
Important: If change the above expression to
we will get an XML parsing error. Why? This i because XML parser will think
that <1
is the beginning of the new XML element. To fix this you could
use two approaches:
1.Present your expression as CDATA
<AdditionalTerm>
<![CDATA[
CellType<1 ? 0.01*F : 0.15*F
]]>
</AdditionalTerm>
In this case XML parser will correctly interpret the expression enclosed
between <![CDATA[
and ]]>
.
2. Replace XML using equivalent Python syntax - see (http://pythonscriptingmanual.readthedocs.io/en/latest/replacing_cc3dml_with_equivalent_python_syntax.html) in which case you would code the above XML element as the following Python statement:
DiffusionDataElmnt\_2.ElementCC3D('AdditionalTerm', {}, 'CellType<1 ? 0.01*F : 0.15*F')
The moral from this story is that if like to use muParser in the XML file make sure to use this general syntax:
<AdditionalTerm>
<![CDATA[
YOUR EXPRESSION
]]>
</AdditionalTerm>
One thing to remember is that computing time of the additional term depends on the level of complexity of this term. Thus, you might get some performance degradation for very complex expressions coded in muParser
Similarly as in the case of FlexibleDiffusionSolverFE
we may use
<AutoscaleDiffusion>
tag tells CC3D to automatically rescale diffusion
constant. See section FlexibleDiffusionSolver
or the Appendix
for more
information.
Steady State diffusion solver¶
Often in the multi-scale simulations we have to deal with chemicals
which have drastically different diffusion constants. For slow diffusion
fields we can use standard explicit solvers (e.g. FlexibleDiffusionSolverFE
)
but once the diffusion constant becomes large
the number of extra calls to explicit solvers becomes so large that
solving diffusion equation using Forward-Euler based solvers is simply
impractical. In situations where the diffusion constant is so large that
the solution of the diffusion equation is not that much different from
the asymptotic solution (i.e. at \(t=\infty\)) it is often more convenient to use
SteadyStateDiffusion
solver which solves Helmholtz equation:
where \(F\) is a source function of the coordinates - it is an input to the equation, \(k\) is decay constant and \(c\) is the concentration. The \(F\) function in CC3D is either given implicitly by specifying cellular secretion or explicitly by specifying concentration \(c\) before solving Helmholtz equation.
The CC3D stead state diffusion solvers are stable and allow solutions for large values of diffusion constants. The example syntax for the steady-state solver is shown below:
<Steppable Type="SteadyStateDiffusionSolver2D">
<DiffusionField Name="INIT">
<DiffusionData>
<FieldName>INIT</FieldName>
<DiffusionConstant>1.0</DiffusionConstant>
<DecayConstant>0.01</DecayConstant>
</DiffusionData>
<SecretionData>
<Secretion Type="Body1">1.0</Secretion>
</SecretionData>
<BoundaryConditions>
<Plane Axis="X">
<ConstantValue PlanePosition="Min" Value="10.0"/>
<ConstantValue PlanePosition="Max" Value="5.0"/>
</Plane>
<Plane Axis="Y">
<ConstantDerivaive PlanePosition="Min" Value="0.0"/>
<ConstantDerivaive PlanePosition="Max" Value="0.0"/>
</Plane>
</BoundaryConditions>
</DiffusionField>
</Steppable>
The syntax is is similar (actually, almost identical) to the syntax of
the FlexibleDiffusionSolverFE
. The only difference is that while
FlexibleDiffusionSolverFE
works in in both 2D and 3D users need to
specify the dimensionality of the steady state solver. We use
<Steppable Type="SteadyStateDiffusionSolver2D">
for 2D simulations when all the cells lie in the xy
plane and
<Steppable Type="SteadyStateDiffusionSolver">
for simulations in 3D.
Note
We can use Python to control secretion in the steady state solvers but
it requires a little bit of low level coding. Implementing secretion
in steady state diffusion solver is different from “regular” Forward
Euler solvers. Steady state solver takes secretion rate that is
specified at t=0
and returns the solution at t=∞
. For a large diffusion
constants we approximate solution to the PDE during one MCS by using
solution at`` t=∞``. However, this means that if at each MCS secretion
changes we have to do three things 1) zero entire field, 2) set
secretion rate 3) solve steady state solver. The reason we need to zero
entire field is because any value left in the field at mcs=N
is
interpreted by the solver as a secretion constant at this location at
mcs=N+1
. Moreover, the the secretion constant needs to have negative
value if we want to secrete positive amount of substance - this weird
requirements comes from the fact that we re using 3:sup:`rd` party
solver which inverts signs of the secretion constants.
An example below demonstrates how we control secretion of the steady
state in Python. First we need to include tag <ManageSecretionInPython/>
in the XML definition of the solver:
<Steppable Type="SteadyStateDiffusionSolver2D">
<DiffusionField>
<ManageSecretionInPython/>
<DiffusionData>
<FieldName>FGF</FieldName>
<DiffusionConstant>1.00</DiffusionConstant>
<DecayConstant>0.00001</DecayConstant>
</DiffusionData>
</DiffusionField>
</Steppable>
In Python the code to control the secretion involves iteration over
every pixel and adjusting concentration (which as we mentioned will be
interpreted by the solver as a secretion constant) and we have to make
sure that we inherit from SecretionBasePy
not SteppableBasePy
to ensure
proper ordering of calls to Python module and the C++ diffusion solver.
Note
Make sure you inherit from SecretionBasePy
when you try
to manage secretion in the steady state solver using Python. This will
ensure proper ordering of calls to steppable and to C++ solver code.
Note
Once you use <ManageSecretionInPython/>
tag in the XML
all secretion tags in the SecretionData
will be ignored. In other words,
for this solver you cannot mix secretion specification in Python and
secretion specification in the XML.
def __init__(self, _simulator, _frequency=1):
SecretionBasePy.__init__(self, _simulator, _frequency)
def start(self):
self.field = CompuCell.getConcentrationField \
(self.simulator, "FGF")
secrConst = 10
for x, y, z in self.everyPixel(1, 1, 1):
cell = self.cellField[x, y, z]
if cell and cell.type == 1:
self.field[x, y, z] = -secrConst
else:
self.field[x, y, z] = 0.0
def step(self, mcs):
secrConst = mcs
for x, y, z in self.everyPixel(1, 1, 1):
cell = self.cellField[x, y, z]
if cell and cell.type == 1:
self.field[x, y, z] = -secrConst
else:
self.field[x, y, z] = 0.0
Warning
Notice that all the pixels that do not secrete have to be 0.0 as
mentioned above. If you don’t initialize field values in the
non-secreting pixels to 0.0 you will get wrong results. The above
code, with comments, is available in our Demo suite (Demos/SteppableDemos/SecretionSteadyState
or
Demos/SteppableDemos/SteadyStateDiffusionSolver
).
Fluctuation Compensator Solver Add-On¶
Fluctuation Compensator is an optional PDE solver add-on introduced in v4.1.2 to account for Metropolis surface fluctuations in field solutions in both a feasible and sensible way. The algorithm is based on that which is described by Marée et. al [1] and imposes total mass conservation in all cellular domains and the medium over the spin flips of a Monte Carlo step. The algorithm does not impose advection by domain deformation, but rather homogeneously applies a correction factor in each subdomain to all field solutions of a solver such that the total amount of each simulated species in each subdomain is unchanged over all spin flips.
The Fluctuation Compensator algorithm consists of the following three rules.
- Rule 1 (mass conservation): In each subdomain, the total amount of each species is unchanged over an arbitrary number of spin flips.
- Rule 2 (uniform correction): Corrections applied to values in each subdomain so to impose mass conservation are uniformly applied.
- Rule 3 (Neumann condition): The amount of species in the copying site of a spin flip are exactly copied to the site of the flip.
Note
Fluctuation Compensator is supported in DiffusionSolverFE and ReactionDiffusionSolverFE.
Exactly one Fluctuation Compensator can be attached to each supported solver instance.
An attached Fluctuation Compensator performs corrections on all fields of the solver
to which it is attached. A Fluctuation Compensator can be attached to a solver in
CC3DML using the tag <FluctuationCompensator/>
at the level of field specification.
For example, to attach a Fluctuation Compensator to a simple use-case with
DiffusionSolverFE might look like the following.
<Steppable Type="DiffusionSolverFE">
<FluctuationCompensator/>
<DiffusionField Name="ATTR"/>
</Steppable>
Demos showing basic usage and comparison are available in
Demos/SteppableDemos/FluctuationCompensator
.
Note
PDE solution field values can be modified outside of the solver routines without
invalidating the correction factors of Fluctuation Compensators so long as
they are notified that field values have been modified. The CompuCell3D library has
a convenience method to do exactly this: updateFluctuationCompensators
. Call this
method after modifying field values and before the next PDE solution step to refresh
Fluctuation Compensators according to your changes.
[1] | Marée, Athanasius FM, Verônica A. Grieneisen, and Leah Edelstein-Keshet. “How cells integrate complex stimuli: the effect of feedback from phosphoinositides and cell shape on cell polarization and motility.” PLoS Computational Biology 8.3 (2012). |