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

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:

\begin{eqnarray} \frac{\partial c}{\partial t} = D \nabla^2c+kc+\text{secretion} \end{eqnarray}

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:

\begin{align*} \frac{\partial c}{\partial t} = D \nabla^2c+kc+\text{secretion} + m_dcd + m_fcf \\ \frac{\partial d}{\partial t} = D \nabla^2d+kd+\text{secretion} + n_cdc + n_fdf \\ \frac{\partial f}{\partial t} = D \nabla^2f+kf+\text{secretion} + p_cfc + p_dfd \end{align*}

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:

  1. Set D, Δt, Δx according to your model
  2. If
\begin{eqnarray} D\frac{\Delta t}{\Delta x^2}>0.16 {\text in 3D} \end{eqnarray}

you will need to to call solver multiple times per MCS.

  1. Set <ExtraTimesPerMCS> to N-1 where:
\begin{eqnarray} ND' = D \end{eqnarray}

and

\begin{eqnarray} D\frac{\Delta t/N}{\Delta x^2} < 0.16 {\text in 3D} \end{eqnarray}

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

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:

\begin{eqnarray} \frac{\partial c}{\partial t} = D \nabla^2c+kc+\vec{\nu} \cdot \vec{\nabla} c + \text{secretion} \end{eqnarray}

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:

\begin{eqnarray} \frac{\partial c}{\partial t} = D \nabla^2c+kc+\text{secretion} \end{eqnarray}

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:

\begin{align*} \frac{\partial c_1}{\partial t} = D \nabla^2c_1+kc_1+\text{secretion} + f_1(c_1,c_2,...,c_N, W) \\ \frac{\partial c_2}{\partial t} = D \nabla^2c_2+kc_2+\text{secretion} + f_2(c_1,c_2,...,c_N,W) \\ {\text ...} \\ \frac{\partial c_N}{\partial t} = D \nabla^2c_N+kC_N+\text{secretion} + f_N(c_1,c_2,...,c_N, W) \end{align*}

where W denotes cell type

Let’s consider a simple example of such system:

\begin{align*} \frac{\partial F}{\partial t} = 0.1 \nabla^2F - 0.1H \\ \frac{\partial H}{\partial t} = 0.0 \nabla^2H + 0.1F \end{align*}

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:

f_F = \begin{cases} 0.1F && \text{if CellType=1}\\ 0.51F && \text{otherwise} \end{cases}

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:

\begin{eqnarray} \nabla^2c-kc=F \end{eqnarray}

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).