Approximate simulation of cortical microtubule models using dynamical graph grammars

Dynamical graph grammars (DGGs) are capable of modeling and simulating the dynamics of the cortical microtubule array (CMA) in plant cells by using an exact simulation algorithm derived from a master equation; however, the exact method is slow for large systems. We present preliminary work on an approximate simulation algorithm that is compatible with the DGG formalism. The approximate simulation algorithm uses a spatial decomposition of the domain at the level of the system’s time-evolution operator, to gain efficiency at the cost of some reactions firing out of order, which may introduce errors. The decomposition is more coarsely partitioned by effective dimension (d=0 to 2 or 0 to 3), to promote exact parallelism between different subdomains within a dimension, where most computing will happen, and to confine errors to the interactions between adjacent subdomains of different effective dimensions. To demonstrate these principles we implement a prototype simulator, and run three simple experiments using a DGG for testing the viability of simulating the CMA. We find evidence indicating the initial formulation of the approximate algorithm is substantially faster than the exact algorithm, and one experiment leads to network formation in the long-time behavior, whereas another leads to a long-time behavior of local alignment.


Minimal CMA DGG
The cortical microtubule array (CMA) grammar presented is a minimal set of rules capable of leading to the long-time behavior of network formation.The grammar consists of six rules: two for growth, two for retraction, a collision rule with three different outcomes, and a reversible state change of growing to retracting.We introduce and briefly discuss these in the following subsections.

Positive MT Growth
The positive growth rule models polymerization and effectively elongates the MT.The rule in equation 1 describes MT growth (polymerization), and is a continuous rule.The position of x 2 is updated by solving an ordinary differential equation limited by the dividing length in the direction of the growing node.The function ρgrow (table 4) can be any function of tubulin concentration, but is set to be a linear function with a constant rate of growth v plus and [Y g ] both of which are set in table 3.

Positive MT Overgrowth
Equation 2 is a discrete rule.When a MT segment becomes long enough, a new node is added in between the previous two nodes and splits the edge.This process is called "overgrowth" with respect to the CMA DGG.A sigmoid function provides a fully activated rate when a dividing threshold has been reached.The new node, 3 , is placed at a point controlled by the γ parameter on the line between x 1 and x 2 .An additional constant c can be used to further scale the rate, but for these simulations it is set to unity. ( (2)

Negative MT Retraction
Negative MT retraction is a continuous rule used to model depolymerization.Equation 3 is similar to the growth rule in 1, but instead functions to limit how much the MT segment can retract.The function ρretract (table 4) can be any function of tubulin concentration, but is set to be a linear function with a constant rate of growth v retract and [Y g ] both of which are set in table 3. The limit, L min , is primarily added to prevent instability from round off errors in the ODE solver caused by two points being too close to one another. (3)

Negative MT Undergrowth
The rule in equation 4 is a discrete rule used to shorten a MT.In the context of the CMA DGG, the process of shortening is "undergrowth".The rule in equation 2 inserts more intermediate nodes into the system and the negative MT undergrowth rule removes them.
Here, a Heaviside function H scaled by a constant c set to unity is used.When the retraction node gets close enough to its intermediate neighbor, the rate function activates.If this rule is selected to occur, the neighboring intermediate node is removed.

MT Collision
The collision rule is more combinatorially complex than the rest of the rules, having three alternative discrete outcomes.The LHS of the grammar rule can be seen in equation 5.It represents the case where a growing end of a MT is close enough to the interior of another MT.
Each of the three discrete outcomes makes use of locally defined rule specific parameters from the where clause in equation 6. Normally, the where clause would follow the with clause like in equation 8, but it is separated out for clarity.In equation 6, I 1 , I 2 are intersection points of lines parameterized by α 1 , α 2 ∈ (−∞, ∞) starting at the point x 2 in the directions u 1 , u 2 .The parameter α 1 indicates where the growing edge ( 1 , 5 ) intersects edge ( 1 , 2 ) and α 2 indicates where the growing edge ( 1 , 5 ) intersects edge ( 2 , 3 ). where The functions β 1 , β 2 are used to calculate the minimum distance to intersection along the respective segments.Since the growing edge ( 1 , 5 ) could intersect ( 1 , 2 ) or ( 2 , 3 ), we must find the closest of the two.Hence, we use β = min(β 1 , β 2 ).
The first discrete outcome is zippering (equation 7), with additional parameters defined in the where clause in equation 6.The rate function primarily depends on the incoming MT being at a critical angle and the incoming MT being on a collision trajectory with another MT in the reaction radius.If these conditions are met, the rate function becomes non-zero and increases as the collision gap narrows.The rewrite on the right hand side (RHS) attaches the growing end to the segment and transforms it into a zippering node.
The second case is junction formation (equation 8), with additional parameters defined in the where clause in equation 6.The rate function, like zippering, depends on the critical angle and the collision trajectory.If the conditions are met, the propensity increases and a junction is formed if the rule is selected to occur.The rewrite effectively crosses the growing MT over the intermediate MT segment. where The third discrete outcome is catastrophe (equation 9), with additional parameters defined in the where clause in equation 6.In this case, we have the propensity increasing as long as the MT is on a collision course with the nearby MT segment.The rewrite causes the incoming MT to change from a growing to a retracting state, simulating catastrophe.

State Changes
The state change rule in equation 10 is a bi-directional discrete rule.Growing ends become retracting ends at a rate of ρretract←growth .Retracting ends become growing ends at a rate of ρgrowth←retract .These rates are typically chosen to be near zero, so that they occur less frequently than the other discrete rules [1].More rates can be found in [2]. ( 2 Calculation of correlation length difference z-score The first and second experiments produced angle correlation lengths as detailed in the main text of ξ 1 ±σ ' = 1.34±0.039and ξ 2 ±σ 2 = 3.14±0.06.We estimate the z-score for the difference This is larger than the Z = 5 standard deviations sometimes used for high certainty because the corresponding p-value for a Gaussian distribution is roughly 3 × 10 −7 .Such p-values are ≪ 0.05, a commonly used threshold for "statistical significance".Of course, we may not have a Gaussian distribution.The worst case distributions however must still obey the Chebyshev inequality under which the p-value must be at most 1 Z 2 ≃ .003≪ .05 or even 1 Ẑ2 ≃ .015≪ .05, in an extremely conservative calculation.

Table 1 :
Graph Node Type Symbol TableTable1 lists the symbols used to represent nodes types in the DGG CMA rules.Table2lists the commonly used rule parameters in the grammar, and includes a brief description of their meaning.In table 3 we do the same, but for model parameters.Finally, table 4 lists commonly used functions and a description as well.

Table 2 :
Rule Parameter Definitions

Table 3 :
Model Parameter Definitions

Table 4 :
Function Descriptions The two experiments resulted in different numbers of microtubule segments.If we zoom in on the first experiment window by a factor of 1.6 this population becomes equal, although other parameters of the model become different.So a very generous interpretation of ξ 1 is to multiply it by 1.6 yielding ξ1 ± σ1 = 1.34 * 1.6 ± 0.039 * 1.6 and a z-score of