Chemotaxis

This tutorial will guide you through the steps of adding a chemotactic force to the cells in your simulation. Requirements for chemotaxis include a population of agents, and a chemical field of a chemoattractant. The agents will move along the chemical gradient toward the chemoattractant.

In this example, a simple method is used to calculate the chemotactic force $\triangle F$ on a particle $p$ which responds to a chemoattractant $C$, as follows:

1. Create a random unit vector ${\vec {c}}$, as a potential chemotactic force acting on $p$.
2. In the direction of ${\vec {c}}$, sample $C$ ahead of $p$, as $C^{+}$, and behind $p$, as $C^{-}$.
3. Calculate the chemotactic force using: $\triangle F=\lambda (C^{+}-C^{-})$
4. If $\triangle F>0$, then apply the force as: $\triangle F\cdot {\vec {c}}$

Here, $\lambda$ is a constant scaling parameter that represents the strength of the chemotactic force along the chemical gradient. In code, this parameter will be defined as a global constant value.

disp = mechIntrctData.force;

/* forward and backward direction vectors */
VReal fwdDir;
VReal bckDir;
REAL mag = 0.0;

/* calculate a random forward vector, and its magnitude */
for (S32 dim = 0; dim < DIMENSION; dim++) {
  fwdDir[dim] = -0.5 + Util::getModelRand(MODEL_RNG_UNIFORM);
  mag += fwdDir[dim] * fwdDir[dim];
}
mag = sqrt(mag);

/* normalize the forward vector, and create the backward vector */
for (S32 dim = 0; dim < DIMENSION; dim++) {
  fwdDir[dim] /= mag; // normalize
  bckDir[dim] = -fwdDir[dim];
}

/* find the nearest neighbor interface in each direction (see separate function) */
S32 fwdInt = findNearestInterface(vOffset, fwdDir);
S32 bckInt = findNearestInterface(vOffset, bckDir);

/* calculate the offset (-1,0,1) for each neighbor interface */
VIdx fwdOffset;
VIdx bckOffset;
for (S32 dim = 0; dim < DIMENSION; dim++) {
  fwdOffset[dim] = dim == fwdInt ? (fwdDir[fwdInt] > 0 ? 1 : -1) : 0;
  bckOffset[dim] = dim == bckInt ? (bckDir[bckInt] > 0 ? 1 : -1) : 0;
}

/* read the chemoattractant values from each of the two neighbors */
NbrBox<REAL> nbr = v_gridPhiNbrBox[DIFFUSIBLE_ELEM_CHEMOATTRACTANT];
REAL fwdVal = 0.0;
REAL bckVal = 0.0;
if (nbr.getValidFlag(fwdOffset)) {
  fwdVal = nbr.getVal(fwdOffset);
}
if (nbr.getValidFlag(bckOffset)) {
  bckVal = nbr.getVal(bckOffset);
}

/* calculate the chemotactic force, and apply if positive */
S32 agentType = state.getType();
REAL chemForce = A_CELL_CHEMOTAXIS_FORCE_STRENGTH[agentType] * (fwdVal - bckVal);
if (chemForce > 0) {
  for (S32 dim = 0; dim < DIMENSION; dim++) {
    disp[dim] += fwdDir[dim] * chemForce;
  }
}

/* limit the maximum displacement within a single time step */
for( S32 dim = 0 ; dim < DIMENSION ; dim++ ) {
  if( disp[dim] > IF_GRID_SPACING ) {
    disp[dim] = IF_GRID_SPACING;
  }
  else if( disp[dim] < ( IF_GRID_SPACING * -1.0 ) ) {
    disp[dim] = IF_GRID_SPACING * -1.0;
  }
}