diff --git a/docs/models/model2/fungi.tex b/docs/models/model2/fungi.tex index 77e9693..aea8c76 100755 --- a/docs/models/model2/fungi.tex +++ b/docs/models/model2/fungi.tex @@ -40,10 +40,12 @@ \section{Model} This document describes a fungal growth simulation implemented with the -BMX framework. Individual hyphae are modeled as connected cylinders that +BMX framework, which is in turn implmented using the AMReX framework. +Withing BMX, individual fungal hyphae are modeled as connected cylinders that are designed to mimic the behavior of a fungal network that is growing by elongation of tip elements and creation of new segments at the growth tips. -The hyphae cylinders are characterized by a length and a radius. Hyphae +The hyphae cylinders are characterized by a length and a radius, as well as +their connectivity to neighboring hyphae. Hyphae that have reached a maximum size stop growing if they are in an interior location while hyphae representing a growth tip can generate one or two new hyphae at the free growth end of the cylinder. The tip splitting events that @@ -59,8 +61,8 @@ \section{Model} bonding of the hyphae to other hyphae in the fungal network. These flags are discussed in more detail below. Because the hyphae are cylinders that are rotationally symmetric about the cylinder axis, only two angles are needed instead of the -customary three angles used to specify the orientation of an arbitrary shape -in three dimensions. The centers of the cylinder end-caps are located at positions +customary three angles used to specify the orientation of an arbitrary rigid shape +in three dimensions. The centers of the cylinder endcaps are located at positions $\rvec_1$ and $\rvec_2$ and the angles $\phi$ and $\theta$ describe the unit vector $\nhat = (\cos(\phi)\sin(\theta),\sin(\phi)\sin(\theta),\cos(\theta))$ pointing from $\rvec_1$ to $\rvec_2$. For a cylinder of length $h$, the vectors $\rvec$, @@ -70,6 +72,14 @@ \section{Model} \rvec_2 & = & \rvec + \frac{h}{2}\nhat \end{eqnarray*} +The simulation describes the growth of a fungal network on the surface of an agar +growth medium that is saturated with a nutrient solution. The surface of the +agar is parallel to the $xy$ plane and the normal to the $z$-axis. Growth is +generally along the agar surface, although the network can potentially push up +in the $z$-direction. Nutrients and waste products can diffuse in the agar +medium and can enter end exit individual hyphae based on mass transport across +the hyhpae membrane. Material can also exchange between attached hyphae. + \section{Chemistry Model} The chemistry for this fungi simulation is similar to one constructed for looking at bacteria or yeast and consists of a set @@ -78,7 +88,7 @@ \section{Chemistry Model} constituents are considered, $A$, $B$ and $C$. The cell volume is also coupled to these constituents. $A$ and $C$ exist both inside and outside the cell (in the -agar medium). The concentrations outside the cell are denoted with the +agar growth medium). The concentrations outside the cell are denoted with the subscript $out$ while those inside the cell are labeled with the subscript $in$. The compound $B$ is only located inside the cell and is impermeable to the cell membrane. @@ -88,12 +98,12 @@ \section{Chemistry Model} A_{in} &\stackrel{k_2,k_{-2}}{\longleftrightarrow}& B_{in} + C_{in} \\ C_{in} &\stackrel{k_3,k_{-3}}{\longleftrightarrow}& C_{out} \end{eqnarray*} -$A$ is only capable of undergoing a reaction inside the cell. To distinguish +$A$, $B$, and $C$ are only capable of undergoing a reaction inside the cell. To distinguish between concentrations of these components and absolute amounts of component, we use brackets [ ]. In this notation, $[A_{out}]$ is the concentration of $A$ in a grid cell and -$A_{out}$ is the total amount (mass) of $A$ in a grid cell. There is also an -additional equation that converts $B_{in}$ into cell volume. This can be -expressed by the equation +$A_{out}$ is the total amount (mass) of $A$ in a grid cell. + +An additional equation converts $B_{in}$ into cell volume \[ B_{in} \stackrel{k_g,k_v}{\longleftrightarrow} V_{cell} \] @@ -124,7 +134,7 @@ \section{Chemistry Model} k_{3} \end{eqnarray*} These equation are for component $C$. The variable $A_{cell}$ is the surface area of the hyphae, -$V_{cell}$ is the corresponding hyphe volume and $V_{grid}$ is the fluid volume in the grid cell. +$V_{cell}$ is the corresponding hyphae volume and $V_{grid}$ is the fluid volume in the grid cell. Similar equations apply to component $A$. Using these expressions for the change in concentration due to transport across the cell membrane, the @@ -145,8 +155,7 @@ \section{Chemistry Model} These reactions also act as sources and sinks for $A_{out}$ and $C_{out}$ in the diffusion equations that govern transport in the agar growth medium. -An additional equation couples the concentration of $B_{in}$ to cell growth. -The growth is governed by the equation +Hyphae growth depends on the internal concentration of $[B_{in}]$ and is governed by the equation \[ \frac{dV_{cell}}{dt} = V_{cell}k_v[B_{in}] \] @@ -174,6 +183,9 @@ \section{Chemistry Model} goes negative. This can be corrected by decreasing the change in $A$ to a value $\Delta A'$ so that $A$ goes to zero instead. However, this means that the changes in $B$ and $C$ must also be modified to reflect the smaller change in $A$. +It is necessary to check both the the final value of the concentration inside +the hyphae, as well as the concentration in the grid cell holding the hyphae, +both remain positive, so this check must be performed twice for each component. Besides transport in and out of the cell coupled to diffusion in the external growth medium, constituants can also be transferred between hyphae that are directly @@ -362,7 +374,7 @@ \section{Cell forces} (\rvec_{i1}-\rvec_{j1})\cdot\rhvec_j & = & - \tau_i\rhvec_i\cdot\rhvec_j + \tau_j\left|\rhvec_j\right|^2 \end{eqnarray*} -This can easily be solved in most cases for $\tau_i$ and $\tau_j$. +These can easily be solved in most cases for $\tau_i$ and $\tau_j$. After a solution is found, it is necessary to check that the $\tau_i$ and $\tau_j$ satisfy the constraints. If one variable violates the constraints then the solution can be found by setting the values of $\tau_i$ or $\tau_j$ to the limits 0 or 1 and @@ -403,10 +415,23 @@ \section{Cell forces} \] The expressions for the forces on the two sites of segment $j$ are similar. +One additional force is required in order to maintain bonding order. If two +hyphae $i$ and $j$ are bonded at sites $\alpha$ and $\beta$, then a linear +restoring force is used to keep the sites in close proximity to each other +\[ +\fvec_{i\alpha j\beta} = -H(\rvec_{i\alpha}-\rvec_{j\beta}) +\] +The constant $H$ governs the magnitude of the restoring force. This force is +designed to maintain the bonding relationships between the hyphae. +Theoretically, the larger the force, the smaller the separation between bonded +hyphae. However, a very large force requires a correspondingly small timestep to +maintain numerical stability, so the magnitude of $H$ is somewhat restricted by +numerical considerations. + \section{Hydrodynamic Drag} -To model hydrodynamic drag on the hyphae, each of the cites at the ends of the -hyphae is the center of a sphere. Both spheres are the same size -are chosen to represent +To model hydrodynamic drag on the hyphae, each of the sites at the ends of the +hyphae is assumed to be the center of a sphere. Both spheres are the same size +and are chosen to represent the hydrodynamic drag on the segment as forces from other segments push it around. Each sphere is subject to a force, $\fvec_1$, $\fvec_2$, resulting in a translational and rotational motion of the segment. For freely interacting @@ -421,7 +446,7 @@ \section{Hydrodynamic Drag} The $3\times 3$ matrices $\Dmat_{ij}$ are diffusion matrices for hydrodynamically interacting spheres. The details of these matrices can be found in Rotne and Prager\cite{RP}. The vectors $\vvec_1$ and $\vvec_2$ are the -velocities of particles 1 and 2 that represent the filament segment. +velocities of particles 1 and 2 that define the hyphae segment. %For our purposes, it is more convenient to express the equations of motion of the %system in terms of net force and torque on the particle pair. To simplify this, we @@ -537,11 +562,11 @@ \section{Hydrodynamic Drag} \;\;\;\;\mbox{if}\;\;h<2a \end{eqnarray*} All other matrix elements are zero in this orientation. Note that these matrix -elements are continuous at the contact value of $h=2a$. This formula can be used to +elements are continuous at the contact value of $h=2a$. These formulae can be used to calculate the diffusion tensor in the frame of reference where the cylinder is located along the $x$ axis and then the tensor can be rotated back to the actual orientation. -Two remaining details are to choose is to choose a value for the radius $a$ and +Two remaining details are to choose a value for the radius $a$ and to account for the fixed separation of the ends of the segment as it translates and rotates. For these simulations we choose the radius $a$ so that the total area of the two spheres is equal to the @@ -737,7 +762,7 @@ \section{Maintaining Connectivities Between Segments} A major function of the simulation is maintaining these connectivities and it is complicated to a significant degree by the behavior of the loops over neighbors. Generally, when new segments are created, they are spawned from a parent segment. At the moment of creation, -the permanent copy of both the parent and the newly created child segment are available to +the permanent copy of both the parent and the newly created child segment are available to the program and modifications to their data will be preserved. Thus, when a growth tip generates a new tip or splits into two new tips, all the permanent information for both the parent segment and the new segments is available and all bonding information can be updated @@ -753,7 +778,8 @@ \section{Maintaining Connectivities Between Segments} internal data so that it it is bonded to A and B. However, at this point B's bonding data does not list C as a bonding partner. This must be fixed up later when the loop to calculate forces is performed. However, to -guarantee that segment B's permanent data is modified, the restriction $i