Difference between revisions of "Chemotaxis"
Line 1: | Line 1: | ||
− | This tutorial will guide you through the steps of adding a chemotactic force to the cells in your simulation. | + | This tutorial will guide you through the steps of adding a chemotactic force to the cells in your simulation. Existing requirements for chemotaxis include a population of agents, a chemical field to contain a chemoattractant, and a method of adding/moving chemoattractant concentrations throughout the field (such as secretion and diffusion with decay). The agents will move along the chemical gradient toward the chemoattractant. |
+ | |||
+ | This page focuses on the approach and implementation specific to the chemotactic force, which samples the chemical field for its chemoattractant concentrations. | ||
+ | |||
+ | == Method == | ||
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: | ||
Line 10: | Line 14: | ||
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. | 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. | ||
+ | == Implementation == | ||
+ | |||
+ | The following code implements the method described above. The code is broken into individual snippets in order to illustrate the steps, but the code may be viewed in whole [https://github.com/djholt/biocellion-chemotaxis/blob/master/model_routine_agent.cpp here]. The code is part of the <code>ModelRoutine::adjustSpAgent()</code> method within <tt>model_routine_agent.cpp</tt>. | ||
+ | |||
+ | |||
+ | Calculate a random, normalized forward vector, and a backward vector in the opposite direction. | ||
<nowiki> | <nowiki> | ||
− | |||
− | |||
/* forward and backward direction vectors */ | /* forward and backward direction vectors */ | ||
VReal fwdDir; | VReal fwdDir; | ||
Line 25: | Line 33: | ||
mag = sqrt(mag); | mag = sqrt(mag); | ||
− | /* normalize the forward vector, and | + | /* normalize the forward vector, and set the backward vector */ |
for (S32 dim = 0; dim < DIMENSION; dim++) { | for (S32 dim = 0; dim < DIMENSION; dim++) { | ||
fwdDir[dim] /= mag; // normalize | fwdDir[dim] /= mag; // normalize | ||
bckDir[dim] = -fwdDir[dim]; | bckDir[dim] = -fwdDir[dim]; | ||
− | } | + | }</nowiki> |
+ | |||
− | /* find the nearest neighbor interface in | + | Determine which two neighbor boxes (of the 26 neighbors) intersect with the forward and backward vectors, and then read the forward and backward chemoattractant values. |
+ | <nowiki> | ||
+ | /* find the nearest neighbor interface (x, y, or z) in both the | ||
+ | forward and backward directions (using a separate function) */ | ||
S32 fwdInt = findNearestInterface(vOffset, fwdDir); | S32 fwdInt = findNearestInterface(vOffset, fwdDir); | ||
S32 bckInt = findNearestInterface(vOffset, bckDir); | S32 bckInt = findNearestInterface(vOffset, bckDir); | ||
− | /* calculate the offset (-1,0,1) for each | + | /* calculate the offset ([-1,0,1] in all dimensions) for each of the two neighbors */ |
VIdx fwdOffset; | VIdx fwdOffset; | ||
VIdx bckOffset; | VIdx bckOffset; | ||
Line 52: | Line 64: | ||
if (nbr.getValidFlag(bckOffset)) { | if (nbr.getValidFlag(bckOffset)) { | ||
bckVal = nbr.getVal(bckOffset); | bckVal = nbr.getVal(bckOffset); | ||
− | } | + | }</nowiki> |
+ | |||
− | /* calculate the chemotactic force, and apply if positive */ | + | Calculate the chemotactic force, apply it (if it's positive), and finally limit the force if it moves the particle more than the width of a single box. |
− | + | <nowiki> | |
− | REAL chemForce = A_CELL_CHEMOTAXIS_FORCE_STRENGTH[ | + | /* calculate the chemotactic force, and apply it if positive */ |
+ | disp = mechIntrctData.force; | ||
+ | REAL chemForce = A_CELL_CHEMOTAXIS_FORCE_STRENGTH[state.getType()] * (fwdVal - bckVal); | ||
if (chemForce > 0) { | if (chemForce > 0) { | ||
for (S32 dim = 0; dim < DIMENSION; dim++) { | for (S32 dim = 0; dim < DIMENSION; dim++) { | ||
Line 71: | Line 86: | ||
disp[dim] = IF_GRID_SPACING * -1.0; | disp[dim] = IF_GRID_SPACING * -1.0; | ||
} | } | ||
− | } | + | }</nowiki> |
− | </nowiki> | + | |
+ | == Steps == | ||
+ | |||
+ | # A | ||
+ | # B |
Revision as of 04:05, 16 December 2016
This tutorial will guide you through the steps of adding a chemotactic force to the cells in your simulation. Existing requirements for chemotaxis include a population of agents, a chemical field to contain a chemoattractant, and a method of adding/moving chemoattractant concentrations throughout the field (such as secretion and diffusion with decay). The agents will move along the chemical gradient toward the chemoattractant.
This page focuses on the approach and implementation specific to the chemotactic force, which samples the chemical field for its chemoattractant concentrations.
Method
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.
Implementation
The following code implements the method described above. The code is broken into individual snippets in order to illustrate the steps, but the code may be viewed in whole here. The code is part of the ModelRoutine::adjustSpAgent()
method within model_routine_agent.cpp.
Calculate a random, normalized forward vector, and a backward vector in the opposite direction.
/* 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 set the backward vector */ for (S32 dim = 0; dim < DIMENSION; dim++) { fwdDir[dim] /= mag; // normalize bckDir[dim] = -fwdDir[dim]; }
Determine which two neighbor boxes (of the 26 neighbors) intersect with the forward and backward vectors, and then read the forward and backward chemoattractant values.
/* find the nearest neighbor interface (x, y, or z) in both the forward and backward directions (using a separate function) */ S32 fwdInt = findNearestInterface(vOffset, fwdDir); S32 bckInt = findNearestInterface(vOffset, bckDir); /* calculate the offset ([-1,0,1] in all dimensions) for each of the two neighbors */ 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, apply it (if it's positive), and finally limit the force if it moves the particle more than the width of a single box.
/* calculate the chemotactic force, and apply it if positive */ disp = mechIntrctData.force; REAL chemForce = A_CELL_CHEMOTAXIS_FORCE_STRENGTH[state.getType()] * (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; } }
Steps
- A
- B