MATHEMATICAL MODELING OF EFFECT OF MICROTUBULE-TARGETED AGENTS ON MICROTUBULE DYNAMIC INSTABILITY

Microtubule-targeted agents (MTAs), widely used in chemotherapy, are molecules that are able to block cancer cell migration and division. Their effect on microtubule (MT) dynamic instability is measured by their influence on observable parameters of MT dynamics such as growth speed, time-based catastrophe frequency, time-based rescue frequency, etc. In this paper, we propose a new mathematical model that is able to reproduce MT dynamics with an appropriate estimation of the main observable parameters. Using the experimental data on paclitaxel effect in presence of EB proteins, we fitted parameters of the model from several drug concentrations. It enable us to understand which non-observable model parameters are able to reproduce the effect of MTAs and thus to highlight a new potential mechanism of action associated with MTAs effect in presence of EB protein. Introduction Microtubule-targeted agents (MTAs) are currently widely used in chemotherapeutic regimens to treat different types of cancer including breast, lung, myelomas, lymphomas, and leukemias. They have a cytotoxic effect on dividing cells by affecting the M phase. MTAs are divided into two groups, i.e. polymerizing and depolymerizing agents, depending on whether they increase or decrease MT polymer mass at high concentrations. At low concentrations both groups of MTAs impair MT functions and arrest cells in mitosis with little influence on polymerized mass [JW98]. It has been proposed that the common cytotoxic mechanism of all MTAs is kinetic stabilization (inhibition) of MT dynamics (growing and shortening) [JW98]. The suppression of MT dynamics by MTAs, such as paclitaxel, has been demonstrated both in vitro [DWJ95] and in cells [JW98]. More recently, the paradoxical effects of MTAs on MT dynamics in vitro, in presence of EB protein, has been revealed [PHM12], [MKD13]. Our aim is to better ∗ The program is funded thanks to the support of the AMIDEX PROJECT (No ANR-11-IDEX-0001-02) funded by the “INVESTISSEMENTS D’AVENIR” French Government program, managed by the French National Research Agency (ANR), and the support INSERM PLAN CANCER No PC201418. 1 Aix Marseille Univ, CNRS, Centrale Marseille, I2M, Marseille, France 2 Aix-Marseille Univ, INSERM, CRO2, Marseille, France c © EDP Sciences, SMAI 2018 Article published online by EDP Sciences and available at https://www.esaim-proc.org or https://doi.org/10.1051/proc/201862186201 2 ESAIM: PROCEEDINGS AND SURVEYS understand this action of MTAs, such as paclitaxel, in different experimental conditions (i.e. drug concentration, presence/absence of EBs etc ...) using mathematical modeling. We want to study how the observable effect on ”macroscopic” level can be produced by the MTA on the ”microscopic” one. In particular, we reproduce the action of paclitaxel at various concentrations on MT dynamics in presence of EB protein through the model parameters. We use the improved model based on the continuous model published by P. Hinow et al [HRT09]. In our previous work, we incorporated into the model [HRT09] a new variable a, called the “MT age” of the population of growing MTs. The value of the variable a provides, in fact, information on the MT growth lifetime. It enables us to obtain the reliable estimations of average timeand distance-based catastrophe frequencies in our model [BGH17]. In the present work, we use the same idea to evaluate the time-based rescue frequency. Indeed, we structure MTs in shortening state by the “shortening age” a. We use the information on the shrinkage lifetime to obtain the estimation on time-based rescue frequency. Thus, our new model of MT dynamics capture main measurable indicators frequently provided in literature: MT growth and depolymerization speeds, time-based catastrophe and rescue frequencies. This allows us to reproduce, through the model parameters, the MTA effect on MT dynamics described by a change in these four indicators for various drug concentrations. 1. Continuous model of microtubule dynamics 1.1. Model The model is based on the continuous approach presented in [HRT09]. The time evolution of concentrations of free GTPand GDP-tubulin are described by the functions of time t, p(t) and q(t), respectively. We denote the densities of the populations of growing and shrinking MTs by the functions t → u(t, ·, ·, ·) and t → v(t, ·, ·), respectively. We structure growing MTs by the length x, length of GTP-cap z and MT age a. The variable a has a dimension of time and starts from zero at the time a MT undergoes a transition from growth to shortening (catastrophe) or vice versa (rescue event). Growing MTs pass in shortening state when they loss their GTP caps, i.e. z = 0. The shrinking MTs are structured by the length x and the MT age ã. Note that, thanks to this new variable, we can collect for each MT, the lifetime and the size variation of each of the growing or shortening events, allowing us to compute accurate indicators of the dynamic instabilites. This variable will also be use to reflect a possible aging phenomenon on the hydrolysis. The dynamics of growing MTs is driven by GTP-tubulin addition rate (the growth rate) γpol(p(t)) and by the rate of GTP hydrolysis γhydro(a). The scheme of MT dynamics is presented on Fig. 1.1. The time evolution of MTs of this population is described by the following transport equation ut + ua + divxz(B(t, a)u) = 0, B(t, a) = ( γpol(p(t)) γpol(p(t))− γhydro(a) ) , a > 0, x > z > 0. (1) This model express the fact that the MT growth speed is a function γpol of the number of GTPtubulin with γpol an increasing function of GTP-tubulin concentration defined in the following way γpol(p) = αpol × 1(p>ps) + ( αpol p(t)− pc ps − pc ) × 1(pc<p<ps). (2) The parameter pc stands for the critical concentration of GTP-tubulin required for polymerization. We assume that the growth speed saturates at high GTP-tubulin concentrations. We denote by the ESAIM: PROCEEDINGS AND SURVEYS 3 Hydrolysis speed γh Growth speed γpol Total length x Stabilizing cap of length z Fragmentation Recycling Polymerization Figure 1. Main features of MT dynamics parameter ps the saturating value of GTP-tubulin concentration and αpol is the maximal growth rate. If the quantity of GTP-tubulin p(t) is between values pc and ps, the growth rate depends linearly on p(t). The transport equation (1) also says that the growth speed of the GTP-Tubulin cap is a balance between the polymerisation and hydrolysis velocities : R(t, a) = γpol(p(t))− γhydro(a). (3) In our model, the GTP-hydrolysis mechanism is vectorial, which means that it propagates as a front along a MT body. The speed of the front dividing GTPand GDP-tubulin dimers γhydro(a), which we call GTP-hydrolysis rate is an increasing function of age a γhydro(a) = γ young hydro × 1(ac<a<as) + γ old hydro × 1(as<a). (4) The parameter ac is a critical age that corresponds to the the time, observed in the experiment, necessary for the freshly incorporated GTP-tubulin dimer to be hydrolyzed. We set the GTPhydrolysis to be equal to zero for MT with an age less than ac. For the MT with the age a which is ac < a < as, we set the GTP-hydrolysis rate to be equal to the parameter γ young hydro . We call these MTs “young ones”. We call MTs with the age a > as “old MTs” and set for them higher value of the GTP-hydrolysis rate γ hydro > γ young hydro . The nucleation process appears as a boundary condition on Γ1 = {a > 0, x > 0, z = x} γhydro(a)u(t, a, x, x) = N (p)ψ(x)Θ(a), N (p) = μp, t ∈ [0,∞), where μ is a nucleation rate. We introduce here two non-negative functions Θ(a) and ψ(x) such that ∫ ∞ 0 Θ(a)da = 1 and ∫ ∞ 0 xψ(x)dx = 1 (5) with Θ and ψ supported on supp Θ = (0, a0), and supp ψ = (0, x0), respectively. The function ψ(x) defines the normalized MT length distribution. 4 ESAIM: PROCEEDINGS AND SURVEYS Rescue events appear as a boundary condition on Γ2 = {a > 0, x > 0, z = 0} R(t, a)u(t, a, x, 0) = Θ(a)λ ∞ ∫ 0 1(ã>ares)v(t, ã, x)dã, if R(t,a) > 0, t ∈ [0,∞). (6) Here, λ characterizes the propensity for shrinking MTs to be rescued. The threshold ares is used to restrict the age ã of shrinking MTs that can return to the polymerization state. We finally supply the system with a boundary condition on Γ3 = {(a, x, z) s.t. a = 0, x > 0, z > 0} compatible with other boundary conditions u(t, 0, x, z) = 0, t ∈ [0,∞). The dynamics of the population of shrinking MTs is described by the transport equation vt(t, ã, x) + divax(Bvv) = Iu→v − Iv→u, Bv = ( 1 −γdepol ) , (7) endowed with v(t, 0, x) = 0, t > 0, x > 0. (8) By Iu→v and Iv→u we denote influx and outflux to the domain of growing MTs that have the form


Introduction
Microtubule-targeted agents (MTAs) are currently widely used in chemotherapeutic regimens to treat different types of cancer including breast, lung, myelomas, lymphomas, and leukemias. They have a cytotoxic effect on dividing cells by affecting the M phase.
MTAs are divided into two groups, i.e. polymerizing and depolymerizing agents, depending on whether they increase or decrease MT polymer mass at high concentrations. At low concentrations both groups of MTAs impair MT functions and arrest cells in mitosis with little influence on polymerized mass [JW98]. It has been proposed that the common cytotoxic mechanism of all MTAs is kinetic stabilization (inhibition) of MT dynamics (growing and shortening) [JW98]. The suppression of MT dynamics by MTAs, such as paclitaxel, has been demonstrated both in vitro [DWJ95] and in cells [JW98]. More recently, the paradoxical effects of MTAs on MT dynamics in vitro, in presence of EB protein, has been revealed [PHM + 12], [MKD + 13]. Our aim is to better understand this action of MTAs, such as paclitaxel, in different experimental conditions (i.e. drug concentration, presence/absence of EBs etc ...) using mathematical modeling. We want to study how the observable effect on "macroscopic" level can be produced by the MTA on the "microscopic" one. In particular, we reproduce the action of paclitaxel at various concentrations on MT dynamics in presence of EB protein through the model parameters.
We use the improved model based on the continuous model published by P. Hinow et al [HRT09]. In our previous work, we incorporated into the model [HRT09] a new variable a, called the "MT age" of the population of growing MTs. The value of the variable a provides, in fact, information on the MT growth lifetime. It enables us to obtain the reliable estimations of average time-and distance-based catastrophe frequencies in our model [BGH + 17].
In the present work, we use the same idea to evaluate the time-based rescue frequency. Indeed, we structure MTs in shortening state by the "shortening age" a. We use the information on the shrinkage lifetime to obtain the estimation on time-based rescue frequency. Thus, our new model of MT dynamics capture main measurable indicators frequently provided in literature: MT growth and depolymerization speeds, time-based catastrophe and rescue frequencies. This allows us to reproduce, through the model parameters, the MTA effect on MT dynamics described by a change in these four indicators for various drug concentrations.

Model
The model is based on the continuous approach presented in [HRT09]. The time evolution of concentrations of free GTP-and GDP-tubulin are described by the functions of time t, p(t) and q(t), respectively. We denote the densities of the populations of growing and shrinking MTs by the functions t → u(t, ·, ·, ·) and t → v(t, ·, ·), respectively. We structure growing MTs by the length x, length of GTP-cap z and MT age a. The variable a has a dimension of time and starts from zero at the time a MT undergoes a transition from growth to shortening (catastrophe) or vice versa (rescue event). Growing MTs pass in shortening state when they loss their GTP caps, i.e. z = 0. The shrinking MTs are structured by the length x and the MT ageã. Note that, thanks to this new variable, we can collect for each MT, the lifetime and the size variation of each of the growing or shortening events, allowing us to compute accurate indicators of the dynamic instabilites. This variable will also be use to reflect a possible aging phenomenon on the hydrolysis.
The dynamics of growing MTs is driven by GTP-tubulin addition rate (the growth rate) γ pol (p(t)) and by the rate of GTP hydrolysis γ hydro (a). The scheme of MT dynamics is presented on Fig. 1.1. The time evolution of MTs of this population is described by the following transport equation This model express the fact that the MT growth speed is a function γ pol of the number of GTPtubulin with γ pol an increasing function of GTP-tubulin concentration defined in the following way The parameter p c stands for the critical concentration of GTP-tubulin required for polymerization. We assume that the growth speed saturates at high GTP-tubulin concentrations. We denote by the parameter p s the saturating value of GTP-tubulin concentration and α pol is the maximal growth rate. If the quantity of GTP-tubulin p(t) is between values p c and p s , the growth rate depends linearly on p(t). The transport equation (1) also says that the growth speed of the GTP-Tubulin cap is a balance between the polymerisation and hydrolysis velocities : In our model, the GTP-hydrolysis mechanism is vectorial, which means that it propagates as a front along a MT body. The speed of the front dividing GTP-and GDP-tubulin dimers γ hydro (a), which we call GTP-hydrolysis rate is an increasing function of age a γ hydro (a) = γ young hydro × 1 (ac<a<as) + γ old hydro × 1 (as<a) .
The parameter a c is a critical age that corresponds to the the time, observed in the experiment, necessary for the freshly incorporated GTP-tubulin dimer to be hydrolyzed. We set the GTPhydrolysis to be equal to zero for MT with an age less than a c . For the MT with the age a which is a c < a < a s , we set the GTP-hydrolysis rate to be equal to the parameter γ young hydro . We call these MTs "young ones". We call MTs with the age a > a s "old MTs" and set for them higher value of the GTP-hydrolysis rate γ old hydro > γ young hydro . The nucleation process appears as a boundary condition on where µ is a nucleation rate. We introduce here two non-negative functions Θ(a) and ψ(x) such that with Θ and ψ supported on supp Θ = (0, a 0 ), and supp ψ = (0, x 0 ), respectively. The function ψ(x) defines the normalized MT length distribution.
Rescue events appear as a boundary condition on Γ 2 = {a > 0, x > 0, z = 0} Here, λ characterizes the propensity for shrinking MTs to be rescued. The threshold a res is used to restrict the ageã of shrinking MTs that can return to the polymerization state. We finally supply the system with a boundary condition on Γ 3 = {(a, x, z) s.t. a = 0, x > 0, z > 0} compatible with other boundary conditions The dynamics of the population of shrinking MTs is described by the transport equation endowed with v(t, 0, x) = 0, t > 0, x > 0.
By I u→v and I v→u we denote influx and outflux to the domain of growing MTs that have the form The " − " operator is given in (10).
The time evolution of GTP-and GDP-tubulin concentrations p and q, respectvely, are described by following ordinary differential equations where I p→u =

Macroscopical output values of the model
In the Table 1 you can find the notations of output macroscopical indicators of MT dynamics and their estimations. There are two groups of output values: observable and non-observable in the experiment. Note that in description of distance-based catastrophe frequency, we use the expression

Notation Mathematical expression Meaning
Values observable in experiment The growth (or polymerization) speed (or rate)

Decoration time
Values non-observable in experiment or hard to estimate The number of MTs in growth state The number of MTs in depolymerization (or shrinking or shortening) state Average GTP-hydrolysis rate  Decoration time T deco is the time of existence of EB1 binding sites [PHM + 12]. This formula of T deco is provided by the biological estimation. It is based on the assumption that the value of the growth rate, which is observable parameter, is probably close to the value of GTPhydrolysis rate which is non-observable. In our model, more appropriate alternative formula could be L av cap (t)/γ av hydro . However, we choose the formula using as a biological estimation to be able to compare our results with experimental ones, in the future.

Lemma 1 (Conservation of tubulin)
The total quantity of tubulin in the system, which is a sum of free GTP-and GDP-tubulin, and the quantities of tubulin contained in MTs in both polymerization and depolymerization states,

Model parameters
Our numerical results provided in [B17] showed that we can divide our parameters on two groups. In the first group are the parameters that mainly affect output "kinetic" indicators such as the growth speed, the time-based and distance-based catastrophe frequencies, the time-based and distance-based rescue frequencies without strong influence on proportion of polymerized-to-free tubulin mass. In another group are the parameters that poorly affect the "kinetic" indicators but strongly change the proportion of polymerized-to-free tubulin mass. In our simulation of paclitaxel effect we will use only the first group of parameters, since we will refer to the biological observations [PHM + 12], describing the effect of MTAs at low concentrations, where few variation in proportion of polymer-to-free tubulin mass were observed. The description of the parameters of the first group is below. Parameters linked to depolymerization.
• The parameter γ depol is directly given by the observable shortening rate.
• The parameter a res which is the time-delay before rescue event regulates the value of the time-based rescue frequency in our model. This time-based rescue frequency can also be adjusted by the parameter λ which is the propensity of rescue. But we choose the parameter a res to regulate the time-based rescue frequency, since a res is simpler to interpret from the biological point of view than the parameter λ. Moreover, if we get the information on both shortening rate and the distance-based rescue frequency then we get an upper bound of a res . In the case, when the distance based rescue frequency is not much affected by the drug the increase of the shortening rate probably lead to the decrease of the parameter a res .
Thus, by the γ depol and a res we can fix the observable parameters shortening rate and time-based rescue frequency at the values that we get from the experiment. Parameters linked to the GTP-hydrolysis. The parameters defining the GTP-tubulin hydrolysis rate can not be easily estimated in the experiment. These parameters are the initial time delay of the GTP-tubulin hydrolysis a c , the two GTP-hydrolysis rates γ young hydro and γ old hydro , and the threshold a s defining the age at which MTs undergo the acceleration of GTP-hydrolysis rate from γ young hydro to γ old hydro . Our numerical results reveal that • The parameter a c mainly affects the time-based catastrophe frequency and has a moderate effect on the growth speed. • The parameter of the GTP-hydrolysis rate γ young hydro has more effect on the growth speed than on time-based catastrophe frequency.
• The parameter of the GTP-hydrolysis rate γ old hydro has the similar effect as γ young hydro on the growth speed and on time-based catastrophe frequency in the case (γ * pol > γ young hydro ). In the case (γ * pol < γ young hydro ) this parameter has poor or negligible effect on MT dynamics. • In the case (γ * pol > γ young hydro ) the parameter a s has the similar effect on the growth speed and the time-based catastrophe frequency as the parameter a c . Otherwise, this parameter has a weak or negligible effect on MT dynamics. Both indicators γ * pol and F temp cat change in the same directions if we vary one of the four parameters a c , a s , γ young hydro or γ old hydro . If, two parameters a c and γ young hydro are sufficient to fit the growth rate and the time-based catastrophe frequency, we set γ old hydro = γ young hydro , a s = a c that is equivalent to the one plateau of GTP-hydrolysis rate.
In presence of EB proteins, the experimental data provides information on EB comet length by fluorescence intensity of +end. As EB comet length is close to the GTP-tubulin cap length, therefore, one can estimate the decoration time. In our model, the parameter a c (of GTP-tubulin hydrolysis rate) correlates with this indicator. Thus, in our simulations we can take into account the data on EB fluorescence intensity adjusting the parameter a c .

Meshes
• Let 0 = t 0 < . . . < t n < . . . t N t = T be a discretization of time interval ]0, T [. The nonconstant time step, denoted by dt n = t n+1 − t n is evaluate at each time step to ensure the stability of the scheme. Note that dt n ≤ 0.001 in the numerical results presented in this paper. • Points of discretization in age will evaluate with each time step dt n . On n-th iteration the points of discretization are a n k = a n−1 k−1 + dt n−1 , k = 1, . . . , n. Note that a n k+1 − a n k = dt n−k := ε n k , with a n 0 = 0 and then ε n+1 k+1 = ε n k , ∀k ≥ 1.

Discrete equations
The system (1)-(11) is discretized using an explicit Euler scheme in time, an upwind scheme for the transport terms and order one quadratures for the integral terms insuring the total preservation of tubulin. The time evolution of the quantity of free GTP-tubulin p on a discrete level where The time evolution of the quantity of free GDP-tubulin q on a discrete level where where I n,i u→v = n k=1 ε n k (R n k ) − u n ki1 , endowed with v n 1i = 0, ∀i, and R n k = γ pol (p n ) − γ hydro (a n k ). Discrete equations for u n+1 for squares cells M ij , i = 2, . . . , N x , j = 1, . . . , i − 1, k = 2, . . . , n are given by and for triangular cells M ii , i = 1, . . . , N x , k = 2, . . . , n endowed by the following conditions for all n, i = 1, . . . , N x , j = 1, . . . , i, k = 2, . . . , n u n 1ij = 0, u n ki0 = 0, u n kii+1 = 0, and the boundary conditions γ hydro (a n 1 )u n 1ii+1 = N (p)ψ(x i )Θ n 1 , The boundary condition (6) at the discrete level reads as follows

Biological observations
Consider the data from [PHM + 12] given in Table 2 for the low dose of paclitaxel (1 nM, 10 nM, 100 nM) on MT dynamics in vitro, obtained in condition of presence of 75 nM EB3. We can see that, at all doses of paclitaxel, all observable parameters are changed significantly compared to the control observation. Note that the shortening rate has the wide range of variability. Paclitaxel increased MT dynamic instability by increasing both MT growth and shortening rates. MT growth rates are increased in concentration dependent manner. Paclitaxel increased time-based rescue and catastrophe frequencies, while the space-based catastrophe frequency was poorly affected. The experimental data [PHM + 12] also contains the information on the fluorescence intensity (mean FI in the Table 2) of EB3-GFP comets. To simulate the drug effect we use the information given in Table 2, including the fluorescence intensity of EB-binding site.  Rhomb diagram representation. On Fig. 3 you can see the graph corresponding to the experimental data given in the Table 2. Such kind of a graph, that we call a "rhomb diagram", was firstly presented in [LBD + 14]. It illustrates the action of paclitaxel on four main indicators of MT dynamics: depolymerization rate, growth rate, time-based catastrophe and rescue frequencies. This is a convenient representation to analyze the data on four parameters obtained at different drug concentrations. The diagram allows us visualize the tendency of changing of key parameters, as well as to see if this changing is in correlation with the drug concentration and also to visualize a relative change of the parameter to the control observation. All parameter values are normalized to the values of the control observation. On the Fig. 3 the red rhomb having its' vertices at (−1, 0) (shortening rate), (0, 1) (time-based rescue frequency), (1, 0) (growth speed) and at (0, −1) (timebased catastrophe frequency) corresponds to the control observation. Other rhombs reflect the experimental data obtained in presence of paclitaxel at 1 nM, 10 nM and 100 nM concentrations. Further we will use the rhomb diagrams to represent numerical results of our simulations.

Strategy of simulation of the paclitaxel effect
The aim is to understand which of the model parameters can produce the output values close to ones obtained in the experiment at time T , when the system reaches the steady state.
The successive steps to calibrate the main parameters of the model for each drug concentration are the following: • At first, we fix the parameters linked to depolymerization process.
(1) Fix γ depol according to the value of the shortening rate obtained in the experiment.
(2) Adjust the parameter a res to fit the time-based rescue frequency. • Then fix the parameters of GTP-hydrolysis.
(3) reduce a c in agreement with the fluorescence intensity. This step is due to our observation that the parameter a c is found to regulate the decoration time which, in turn, correlates with fluorescence intensity.
(4) Adjust γ young hydro and if necessary involve the parameters of the second plateau of GTP hydrolysis rate: a s and γ old hydro to fit γ * pol and F temp cat with experimental growth speed and time-based catastrophe frequency. By default, a s = a c and γ old hydro = γ young hydro (i.e. only one plateau is considered), which corresponds to the absence of aging process. This last step might be tricky and requires to perform a big number of numerical tests. This is due to the facts that all parameters a c , a s , γ young hydro and γ old hydro , defining the function of GTPhydrolysis rate, regulate γ * pol and catastrophe frequencies. In particular, γ * pol depends on the average hydrolysis rate which, in turn, depends on interrelation of these four parameters.
Remark. Both γ young hydro and γ old hydro affect more the growth rate than catastrophe frequencies. In contrast, the time delay a c have stronger influence on time-based catastrophe frequency than on growth speed [B17].

Numerical results
Fitting the model parameters. We started from the calibration of model parameters referring to the data corresponding to control observation given in Table 2. The obtained set of parameters can be found in Table 3. In all simulations, we switch off the nucleation process (i.e. we set µ = 0) at time t = 15 min. The initial conditions are the following p(0) = 15 µM, q(0) = 0, u(0, ·, ·, ·) = 0, v(0, ·, ·) = 0.
The simulation performed with the set of parameters given in Table 3 we call "control test". It produces the values of growth speed, time-and distance-based catastrophe frequencies and of timebased rescue frequency which are in strong agreement with the experimental control observations given in Table 2. The output values provided in Table 4 are obtained at time T = 20 min, for which we already observe the steady state.  Table 3. Parameter set reproducing the MT dynamics in vitro of control observation described in [PHM + 12] Observable output indicators Non-observable output indicators  Table 4. Output indicators of MT dynamics obtained for simulated "control observation" corresponding to the parameter set given in Table 3.
Using our strategy we tried to simulate the effect of 1 nM, 10 nM and 100 nM of the paclitaxel described in the Table 2. For this purpose, we change several parameters (including µ = 0) of the control set given in Table 3 at t = 15 min. Thus, in all simulations we examine the influence of the parameter on the same MT population that is formed until T = 15 min.
Numerical results obtained at each step of the strategy and the experimental data, as well, are further represented in rhomb diagrams for each drug concentration. Simulation of 1 nM paclitaxel effect. At the first step, we fix γ depol = 17.9 that has a negligible effect on other parameters (see Fig. 4(a)). After the second step, adjusting the parameter a res , the time-based catastrophe frequency is strongly increased (see Fig. 4(b)). The change of this indicator appeared to be in the same direction as we want to obtain. However, the growth rate slightly reduced and changed in opposite direction to the one provided in the real data. We fix the value of a c at the value which is smaller than in the control test (-10%) similar to the 10% decrease in fluorescence of EB comet [PHM+12]. It results in the slight increase in time-based catastrophe frequency (see Fig. 4(c)). Then by performing a number of tests, we find the following set of parameters: γ depol = 17.9, a res = 0.22, a c = 0.045, a s = 0.11, γ young hydro = 3.5, γ old hydro = 5.   Table 5. Indicators of MT dynamics obtained for simulated effect of 1 nM of paclitaxel.
Simulation of 10 nM paclitaxel effect. After the first step, we observe a weak effect on timebased catastrophe frequency, which is reduced. After the second step, the catastrophe frequency changed in the direction corresponding to the drug effect, however a bit exceeded the required value. The growth speed moderately decreased and became too low compared with the experimental one which is strongly increased (see Fig. 5(d)). We continue to decrease a c in concentration dependent manner (see Fig. 5(c)) that does not influence on the parameters. The found set of parameters is γ depol = 21.3, a res = 0.13, a c = 0.041, γ old hydro = 6.3, a s = 0.

Observable output indicators
Non-observable output indicators  Table 6. Indicators of MT dynamics obtained for simulated effect of 10 nM of paclitaxel.
Simulation of 100 nM paclitaxel effect. After the first step, we observe a strong influence on time-based catastrophe frequency, which is reduced and weak effect on the growth speed which is increased. After the second step of the strategy, we see the opposite effect on these two parameters. The growth speed is slightly reduced and the time-based catastrophe frequency is strongly increased compared to the control one. On the following step, we continue to decrease the parameter a c (see Fig. 6(c)) that increased the time-based catastrophe frequency and do not changed the growth speed. On the fourth steps, we increase γ old hydro and strongly reduce a s . The found set of parameters is γ depol = 28, a res = 0.18, a c = 0.036, a s = 0.04, γ old hydro = 7.  Table 7. Indicators of MT dynamics obtained for simulated effect of 100 nM of paclitaxel.

Observable output indicators Non-observable output indicators
The comparison of the simulated effect for 1 nM, 10 nM and 100 nM and experimental data are given on Fig. 7. The parameter values that were changed from the "control" are given in Table 8 for final simulations of the drug effect.

Conclusions and discussions
We created a new model of MT dynamics allowing to simulate the MTA effect on four measurable indicators of dynamics: a growth rate, a depolymerization rate, a time-based catastrophe frequency and a time-based rescue frequency. We have tested our model by reproducing the paclitaxel effect in presence of EBs. Our model proposes a new concept of the GTP-hydrolysis that is vectorial, age-dependent and delayed. Our analysis showed that the parameter a c as well as parameters linked to MT aging can play a crucial role in regulation of the MT dynamics.
According to the model, we can determine the effect of paclitaxel on non-observable parameters in experiment. We find that in order to fully explain experimental results in presence of EBs paclitaxel might have an effect on the initial time delay of GTP-hydrolysis a c and increase the aging properties of a MT described by parameters a s and γ old hydro in our model. The parameter a s (the time at which the GTP-hydrolysis accelerates) decreases and the second plateau of GTPhydrolysis γ old hydro increases for higher concentrations of paclitaxel in presence of EB proteins. Our results highlight a new potential mechanism of action associated with MTAs effect in presence of EB protein which is induction of MT aging by the GTP hydrolysis during MT growth, including the effect of the initial delay a c . Such effect can be observed on kymographs in EB3-tip tracking assay. Further, we intend to test our strategies on simulation of destabilizing drug effects, such as Vinca Alcaloids.
One of the limits of our model is that it comprises the fixed shortening rate. The fluctuations in shortening rate among one population of MTs can be very strong. In further work, we can represent the depolymerization rate as a fragmentation of MTs to improve our model. We also do not consider in this work the pause state of MTs. It can be important in modeling of drug effects that are known to induce a pause state for a big part of MT population. The pause state can be introduced, for example, by the delay in depolymerization rate.