Difference between revisions of "Chemotaxis"
From Biocellion
Line 3: | Line 3: | ||
In this example, a simple method is used to calculate the chemotactic force <math>\triangle F</math> on a particle <math>p</math> which responds to a chemoattractant <math>C</math>, as follows: | In this example, a simple method is used to calculate the chemotactic force <math>\triangle F</math> on a particle <math>p</math> which responds to a chemoattractant <math>C</math>, as follows: | ||
− | # | + | # Create a random unit vector <math>\vec c</math>, as a potential chemotactic force acting on <math>p</math>. |
− | # | + | # In the direction of <math>\vec c</math>, sample <math>C</math> ahead of <math>p</math>, as <math>C^+</math>, and behind <math>p</math>, as <math>C^-</math>. |
− | # | + | # Calculate the chemotactic force using: <math>\triangle F = \lambda (C^+-C^-)</math> |
+ | # If <math>\triangle F > 0</math>, then apply the force as: <math>\triangle F \cdot \vec c</math> | ||
+ | |||
+ | Here, <math>\lambda</math> 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. | ||
+ | |||
+ | <nowiki> | ||
+ | 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; | ||
+ | } | ||
+ | } | ||
+ | </nowiki> |
Revision as of 23:12, 15 December 2016
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 on a particle which responds to a chemoattractant , as follows:
- Create a random unit vector , as a potential chemotactic force acting on .
- In the direction of , sample ahead of , as , and behind , as .
- Calculate the chemotactic force using:
- If , then apply the force as:
Here, 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; } }