Physical Modeling of Microtubules Network

Microtubules (MT) are highly dynamic tubulin polymers that are involved in many cellular processes such as mitosis, intracellular cell organization and vesicular transport. Nevertheless, the modeling of cytoskeleton and MT dynamics based on physical properties is difficult to achieve. Using the Euler-Bernoulli beam theory, we propose to model the rigidity of microtubules on a physical basis using forces, mass and acceleration. In addition, we link microtubules growth and shrinkage to the presence of molecules (e.g. GTP-tubulin) in the cytosol. The overall model enables linking cytosol to microtubules dynamics in a constant state space thus allowing usage of data assimilation techniques.


Introduction
Microtubules (MTs) are essential for intracellular organization, cell division and cell morphogenesis. It is established that dynamic instability of MTs, with alternating phases of growth and shrinkage separated by catastrophe and rescue events, is regulated by microtubuleassociated proteins (MAPs) [9,4,2,8]. Nevertheless, understanding better the specific molecular interactions at microtubule growing ends recognized by MAPs and deciphering assembly mechanisms in the presence of guanosine triphosphate (GTP) is of major importance in cell biology and for cancer therapy.
Several simulation softwares to visualize the different models of microtubule dynamic instability have been already developed [3,9], but essentially in 2D. In [5], the authors studied MT rigidity using microscopy imaging and physical laws which yielded mechanical properties measures. More recently, Portran et al. showed how MAP65, by modifying MT mechanical properties, regulate the formation of complex MT arrays [8]. In [6], the authors proposed a cytoskeleton simulation model based on derived kinematic laws of mechanics and suitable for fast simulation. However, most of the terms are somewhat empirically designed. Indeed, although this approach is effective for visual simulation, it does not take mass description into account and rather focus on cytoskeleton kinematics than its dynamics. More recently, computation tools that reproduce assembly of bundles made of MTs in the presence of MAP65 have been investigated to mimic in vivo MT bundles [9].
To get a better understanding of how geometrical constraints and interactions with the medium affect the MT organization and dynamics, we propose a new modeling framework which can be combined further with data assimilation techniques. Data assimilation techniques are essentially based on a controlled trade-off between observations of a phenomena and a model accounting for its likely dynamics. Such approach has been for example investigated in [1] for similar but simpler model.
In this paper, we focus on the modeling part of the general data assimilation framework. The model we propose links mechanically the cytoskeleton to the cytosol, and its resulting system allows investigation of an inverse problem. It is composed of two essential parts introduced in the next section: i) the Lagrangian dynamics of MTs (Section 2); ii) the Eulerian dynamics of the cytosol (Section 3).

Microtubule modeling
MTs are discretized by nodes connected by edges and are described either in 2D or 3D. The overall number of nodes N = M P accounts for the size of the state space, while P accounts for the number of nodes for each microtubule, and M the number of microtubules. The set of nodes carries different properties such as positions y ∈ R N d , d ∈ {2, 3}, and velocities u ∈ R N d that can be read for node p i respectively as y i ∈ R d and u i ∈ R d . These properties allows building the following Newtonian equations: where m accounts for the homogeneous mass of nodes, and F f is a force applied to the nodes.
In what follows, we will detail the forces of the model and needed expansions of the system to manage specific behaviors such as rigidity.

Microtubule rigidity
Using the Euler-Bernoulli beam theory, we model the rigidity of MTs on a physical basis. This mechanism relies on a difference between the shape of a given material and its neutral state. In this purpose, we introduce the neutral position state z ∈ R N d accounting for the position of nodes not submitted to any forces.
2.1.1. Elasticity MTs are very rigid material. In [6], MTs are considered as non-elastic by using reshaping algorithm. Although this technique can be justified in practice, it involves an implicit integration potentially expensive from a calculation point of view regarding the number of nodes used for discretization. In order to keep a simple explicit integration, we prefer to consider the MTs as elastic with a high valued spring coefficient. As a counterpart, integration time-step has to be small. By defining c y ij = y j − y i , d y ij = y j − y i and e y ij = c y ij /d y ij the elasticity force reads for each node p i : where s denotes the section of MTs, E is Young modulus and j ↔ i means particle p j is connected to p i . We suppose MTs to be hollow, therefore we can write: s = π(r 2 − r 2 0 ), with r and r 0 being respectively their outer and inner radius.

Bending
Flexibility of MTs is the more tangible mechanism due to rigidity. Before expressing the forces applied on nodes, we detail how angular variation along MTs yields torque. For any node p i connected to nodes p j and p k , the torque is oriented in the direction:  except in the case y i , y j and y k are aligned. We introduce T bending ∈ R N d , the torque produced by angular variations along MTs, reading: where I = π 4 (r 4 − r 4 0 ) is the second moment of area of MTs. Forces produced by torque reads: For d = 3, a torsion component exists, but we ignore it because of the edge modeling of MTs.

Collision
The scale difference between microtubule radius and cytosol domain size ( 10 3 factor) yields numerical difficulties for collision management. A small time step is needed to accurately check whether two microtubules are colliding or not. For simulation ease, the collision is rather modeled as closeness preventing. This kind of modeling yields smooth collision management and easy to tune parameters. The resulting repulsion force is described as follows: where f rep is a reference force and d rep is a reference distance of repulsion. Ideally, for perfect collision management with supposedly hard material, d rep should tend to 0. Instead, we parameter this reference distance accordingly to MT thickness by d rep = αr.

Microtubule growth and shrinkage
In our approach, the number of nodes modeling MTs is constant. The main advantage is that the state space of the model keeps constant too. But as counterpart, managing the extension of curvilinear forms without modifying their shape and with the same discretization can be tedious. We suppose growth/shrinkage velocity to be directed same as microtubule plus-end edge. Considering the neutral position state z of microtubule, as in Fig. 1 (left), intensity of growth velocity is simply projected on nodes with a decreasing ratio such as: where j is the node index inside microtubule M T m , and v growth,m,P is the growth velocity at its plus-end. Such ratio makes the assumption that distance between attached nodes is constant, which is reasonable due to the strong spring coefficient used in the experiments. As a consequence of edge-directed growing velocity, if MTs neutral positions are aligned at t 0 , they will stay so.
In this case, the component γ z jik in (4) is null. Regarding the real microtubule position state y, growth velocity projection applies to u growth almost the same way. But as a difference, the direction of inner nodes velocity varies along microtubule length. For each node, the direction is set as the tangent of the two connected edges, while the same ratio as the neutral state is applied for intensity (see Fig. 1 (right)).  Microtubules index space {j, m} is linked to global index space {i} thanks to the relation: i = h(m, j) = P m + j. That allows us to express growth velocity in the global nodes space such as: v growth,i = v growth,h(m,j) and u growth,i = u growth,h(m,j) . We shall see later how the intensity of growth velocity at plus-end node is determined.

Fluid dragging
MT movements also depend on cytosol dynamics. Various mechanical properties such as pressure or velocity can produce different effects on MT spatial evolution. We only use inner cytosol velocity u(x, t) for viscous dragging of MTs, with x being the position in the cytosol space. Because of fluid-structure interactions, determining accurate drag forces is a tedious task calling for advanced computing fluid dynamics techniques. Fineness of MTs compared to cytosol spatial dimensions and slow velocities yielding a low Reynolds number allows one to make assumption the drag force is simply proportional to fluid velocity. We suppose the force proportional to the dynamic viscosity of cytosol and to the diameter of MTs. The viscous drag force is therefore defined as: F drag,i = 2rηu(y i ).

Cytosol modeling
We propose to link MT growth and shrinkage to the presence of molecules in the cytosol. The growth velocity of MTs depends upon soluble tubulin concentration available for polymerization and GTP-tubulin association and dissociation rates [9,4,2]. The shrinkage velocity is known to be independent of tubulin concentration and is characterized only by the dissociation rate of guanosine diphosphate (GDP) tubulin from the depolymerizing end. In our modeling approach, we do not refer MT growth explicitly to the presence of GTP-tubulin since MT dynamics depend also on numerous properties of the cytosol. In what follows, we only refer to the concentration of a meta-quantity responsible for MT growth and shrinkage. This concentration can be also viewed as a growth potential whose distribution indicates locations of available growth energy. Due to the important role of GTP-tubulin molecules in microtubule polymerization, we call it "pseudo-tubulin" concentration written T (x, t) with the attached unit [T ] = t.m −3 .

"Pseudo-tubulin" modeling
We propose to model evolution of "pseudo-tubulin" using a classical advection-diffusion model with an additional term A accounting for microtubules absorption and release: ∂T with β being the diffusion coefficient of "pseudo-tubulin" and T ref the reference value of T set to 1 t.m −3 . We suppose the absorption term to depend on MT growth velocity at their extremity, as illustrated in Fig. 2. Theoretically the absorption is only localized on this specific position, but due to spatial discretization of the cytosol on a grid, it is necessary to use a kernel convolution, and finally the term reads: where w m values −1 if microtubule m is growing, 1 if shrinking. Also, G σ is the gaussian kernel of standard deviation σ, * is the convolution operator, δ is the Dirac function, and λ stands for the microtubule linear density of "pseudo-tubulin" such as [λ] = t.m −1 . In case "pseudo-tubulin" is considered as α-tubulin and β-tubulin, it is possible to parameter this coefficient with a known value according to their spatial configuration shaping microtubule structure. Otherwise, if T is considered as growing energy, λ can be related to linear energy capacity. Since G σ has to be as close as possible to de Dirac function, σ is set equal to the lowest significant value, being size of grid meshes.

Growth functional
We propose to link microtubules both growth and shrinking occurrence and intensity to "pseudotubulin" concentration T . The intensity of the growth velocity at microtubule plus-end is for instance set as: v growth,m,P = ω λ g(T (y m,P )), with ω standing for the microtubule "pseudo-tubulin" frequency and g(T ) is the intensity function. Their expression changes wether microtubule is growing or shrinking and reads as: ω = ω growth and g We make the assumption that for any growing microtubule m, whenever T (y m,P ) goes below T min , the available growing energy is not sufficient anymore to sustain the growing process. Therefore, the microtubule switch to a catastrophe phase. Oppositely, when microtubule is shrinking, whenever T (y m,P ) goes above T max , the available growing energy is sufficient to start again the growing process.

Final system and experiments
In the end, we obtain a system modeling MT rigidity, growth and interaction with cytosol. The space state describing this system is constant, which is a great advantage for analysis techniques. Finally, the overall system reads:  In our experiments, velocity fields u(x, t) are obtained by simulation of the incompressible Navier-Stokes equations, initialized by null divergence and random curl field. At initialization, "pseudo-tubulin" is randomly distributed around T ref . Hence simulation is able to provide the expected variability. Simulation grid meshes size is set to 1 µm, while integration is performed using fourth order Runge-Kutta scheme with a 2 ms time-step. Results are presented in Figs. 3-4. One can observe our model yields appealing information on cytosol properties. In Fig. 4), simulations have been performed to mimic MT nucleation using MT microseeds adsorbed on a micropatterned glass substrate and observed in TIRF microscopy [7]. This approach is actually used to investigate MT functions using geometrically controlled MT networks in vitro in cell-free assay. Simulation parameters in Table 1 are chosen in agreement with published data.

Conclusion and perspectives
In vivo and in vitro, the MTs form networks resulting from their interplay with other protein partners [4]. The role of mechanical and/or geometrical constraints in the network organization is also important in that concern. Nevertheless, it is difficult to answer simple questions about the mechanics at the level of the whole network. In this paper, we have described a flexible mathematical model for representing MT network dynamics and molecular interactions. This model is able to manage microtubule rigidity, growth and interactions with cytosol. It addresses questions about the mechanisms and the forces developed at the level of MT networks and it has been designed so as to be compliant to analysis techniques based on data assimilation. To go further, we will derive this model in order to apply variational assimilation principles inspired of geophysical analysis. We plan also to combine simulation/data assimilation methods, in vitro experiments to constraint MT dynamics on patterns with specific geometry.