BACKGROUND
fluorescence recovery after photobleaching  was first introduced in  <dig>  <cit>  and is a widely used method to study turnover, transport, diffusion and interaction among biological molecules in living specimens.
the use of frap has been facilitated by the current availability of microscopes equipped with a laser scanning device.
the emergence of fluorescent protein labeling with the green fluorescent protein  and its spectral variants has greatly enhanced frap application.
in a typical frap experiment, a gfp-labeled structure is rapidly and irreversibly photobleached with a high intensity laser, and fluorescence recovery is recorded as a function of time.
the fluorescent molecules then diffuse into the irradiated region, while the non-fluorescent ones diffuse into the unbleached area until equilibrium is reached.
the analysis of fluorescence recovery curves yields the diffusion coefficient and the fraction of free, transiently bound and immobilized molecules.
for a quantitative description of fluorescence recovery dynamics in frap experiments, several theoretical models have been proposed  <cit> .
these include  the pure-diffusion dominant model that considers the recovery rate for weakly bound or free fluorescent molecules, and is defined exclusively by their diffusion;  the effective diffusion model, which describes the recovery kinetics of fluorescent molecules that bind tightly the bleached structure, and is also largely defined by diffusion;  the reaction dominant model, where diffusion is very fast and molecules rapidly equilibrate after the bleach; and  the diffusion phase-binding phase approximation used whenever the contributions of diffusion and binding are coupled.
in another early study, salmon et al.
<cit>  used the frap technique to trace the behavior of fluorescently labeled bovine tubulin injected in sea urchin eggs undergoing the first mitotic division.
the authors provided evidence that the use of exogenous tubulin accurately represents the in vivo situation, and showed that frap data were best fitted by a negative-going exponential function, and that fluorescence recovery had a half-time of approximately 20 s. taking into account that tubulin dimers diffused back into the bleached area within 1 s, the authors concluded that diffusion was not a limiting factor and that the reaction dominant model was correctly interpreting their frap results.
one of the first theoretical description of microtubule growth and degradation  <cit>  was based on data on microtubule dynamics obtained from in vitro experiments  <cit> .
several kinetic models describing microtubule  growth were proposed  <cit> , and one of them, usually referred as the hill’s model  <cit> , is still widely used.
according to this model, mt growth begins with the recruitment of tubulin dimers  at a centrosome-associated nucleation site, followed by the polymerization of additional dimers at a rate constant j
<dig>  or depolymerization at a rate constant j
<dig> .
at any length, the mts may switch from the growth mode to a shrinking mode with a constant rate k
<dig>  or from a shrinking mode to a growth mode with a constant rate k
<dig>  a graphical representation of the hill’s model is reported in fig.
<dig>  this kinetic model was later adapted to describe the polymerization/depolymerization kinetics of microtubules in xenopus egg extracts and to analyze how cyclin a and cyclin b could affect this process  <cit> .
in one of the kinetic regimes , the experimental data were nicely fitting the model, and all of the constants could be defined.
in addition to the bounded state where mts are on average disassembling and j
1
k
2
- j
2
k
<dig> is negative, the authors also described an unbounded state  where mts are on average growing   <cit> .fig.
<dig> kinetic model of mt growth adapted from hill,  <dig>  <cit> .
gi and si are the concentrations of growing and shrinking microtubules containing i number of tubulin dimers.
c <dig> indicates mts with zero length.
j <dig> and j <dig> are the kinetic constants of growth and shrinking, respectively; k <dig>  k <dig> are the constants for rates of growth-to-shrinking and shrinking-to-growth transitions.
growing and shrinking mts are depicted as bars, and tubulin dimers are shown as shaded bars
the hill’s model  was applied to frap-based studies on pre-anaphase b/metaphase spindles of drosophila syncytial embryos expressing gfp-tubulin  <cit> .
the authors concluded that both bounded and unbounded regimes are inadequate to describe the observed frap dynamics.
data analysis indicated that tubulin dimers turn over almost entirely during a single cycle of mt shortening and growth, and consequently the recovery time does not depend on the size or position of the bleached region along the metaphase spindle.
the recovery halftime  after photobleaching is then simply the halftime of this cycle .
this critical analysis served as a starting point for us to find other solutions to the kinetic scheme shown in fig.
<dig>
numerical modeling of a frap experiment  <cit>  shows that the parameters that describe mt turnover can be determined from the modeling data.
however, this approach is based on the mechanical model of brust-mascher et al.
<cit>  and requires equations to be solved for the entire spindle, necessitating the knowledge of multiple mechanical parameters.
as mentioned above, another approach to analyze mt dynamics is based on chemical kinetics description  <cit> .
this approach operates with a smaller number of dynamic parameters and was verified by in vitro microtubule growth experiments.
one of the objectives of the present study is developing a chemical-kinetic description for frap experiments on metaphase spindles.
another objective is the analysis of the dependence of spindle frap on parameters describing microtubule end dynamics.
frap in a steady-state spindle has a simple relationship with the number of mt growing and shrinking ends within the photobleached area.
the dynamics of these ends both with and without specific marks are not easy to evaluate, but can be inferred from the solution of system  for elementary intervals.
thus, our aim is describing the limit transition from the discrete hill’s equations  <cit>  to a continuous equation, and understanding how the kinetic constants of the hill equation compare to those in the continuous equations.
we found a solution for the cauchy problem for partial differential equation   for elementary intervals and used this expression to solve the equation for fluorescence recovery in a steady state spindle.
finally, we applied our model to the experimental data  <cit>  and performed a quantitative analysis of frap recovery time.
RESULTS
limit passage to the continuous model
since we were not able to solve discrete equations describing mt behavior with time, we are using a continuous model.
we denote the mt concentration found in the growth phase as g
i
, where i is the number of polymerized tubulin dimers.
similarly, the concentration of microtubules in the shrinking phase is denoted as s
i
.
then, the kinetic equations for the concentration of tubulin dimer appear as:  <dig> ∂∂tgit=−j1git−gi−1t−k1git+k2sit∂∂tsit=−j2sit−si+1t−k2sit+k1git
we can rewrite equation  as follows, considering Δx is a small coordinate step: ∂∂tgit=−j1ΔxΔxgit−gi−1t−k1git+k2sit∂∂tsit=−j2ΔxΔxsit−si+1t−k2sit+k1git
the values v
1
= j
1
Δx and v
2
= j
2
Δx have a dimension of cm/sec, which correspond to the linear polymerization and depolymerization speeds of a single tubulin dimer.
if Δx is negligible compared to the microtubule size, the above equations can be transformed into pdes, where t is a time and x is a coordinate:  <dig> ∂∂tgxt=−v1∂∂xgxt−k1gxt+k2sxt∂∂tsxt=v2∂∂xsxt−k2sxt+k1gxt
as reported below, this transformation will help us to find a mathematical solution for mt behavior.
stationary regime and bounded state
let us solve the equations for the stationary regime, i.e., when the process does not depend on time, namely when t =  <dig>  in this case, g and s denote the distribution of growing and shrinking ends over the coordinate x, which corresponds to the current length of mts.
v1∂∂xgx=−k1gx+k2sxv2∂∂xsx=−k1gx+k2sx
the discriminant of this system, d=k1v1−k2v <dig> is either positive or zero.
then a characteristic equation: λ2−k2v2−k1v1λ= <dig> has two roots, λ1= <dig> λ2=k2v1−k1v2v1v2=−α
substituting the experimental values of k
1
; k
2
; v
1
; v
<dig> from reference  <cit>  into the expression for λ
<dig> results in a positive α. the general solution of the system  is gx=c1k2v1+c2k2v1e−αxsx=c1k1v1+c2k2v2e−αx where c
<dig> and c
<dig> are arbitrary constants.
if at x = 0; g = g
0
, and upon x tending to infinity g = g
∞, then gx=g∞+g0−g∞e−αxsx=k1k2g∞+g0−g∞v1v2e−αx
using the experimentally observed values of dynamic parameters  <cit>  and arbitrarily assigning g
0
=  <dig> and g
∞
=  <dig>  the system can be visualized as shown in fig.
<dig> fig.
<dig> concentrations of growing and shrinking microtubule ends as a function of coordinate x for a model case
if x is measured in μm, the gk
1
= sk
<dig> ratio holds true for a wide range of large x values, whereas in the case of small x values, which is most pertinent to the experimental situation, an alternative ratio should be applied.
a negative-going exponential function for microtubule length  observed for the discrete form of equations  in the stationary regime  was discussed earlier  <cit>  and was successfully used to interpret the mt length distribution in experiments in vitro.
based on these experimental data  <cit>  and given the finite number of tubulin molecules in the cell, g
∞
>  <dig> is not a viable option and so g
∞ should be set as  <dig>  in this case, a simple ratio gv
1
= sv
<dig> is true and corresponds to a situation where the spindle remains unchanged with time, and the numbers of growing and shrinking microtubules are equal.
analytical solution of the equations for elementary intervals
the system of continuous equations  was first published by dogterom and leibler  <cit> ; the authors stated that it could be solved analytically with the border conditions set at t =  <dig>  and with all mts having zero length.
however, they did not publish the solution.
in the following section we find this solution.
the solution refers to a frap experiment that involves photobleaching of fluorescently labeled mts in a rectangular region of width l perpendicular to the spindle axis, followed by the analysis of the dynamics of fluorescence recovery in the bleached region.
in a real experiment both g  and s  mt ends are distributed uniformly across this region at the initial time-point after bleaching.
let us consider the elementary interval Δ of l, which initially contains g
<dig> and s
<dig> unmarked ends .
the g
<dig> ends will start to produce a marked tubulin track and will move towards the spindle equator from the elementary interval, while the s
<dig> ends will move to the opposite direction.
some g
<dig> ends will then turn into s ends  and some s
<dig> ends will turn into g ends .
to describe the behavior of g
<dig> and s
<dig> ends it is necessary to solve the cauchy problem for .
we thus substituted the x and t variables with α and β, and converted the functions g and s into u and v: α=x−v1tβ=x+v2tgxt=uαβe−αβ+basxt=vαβe−αβ+baa=k1v1+v2b=k2v1+v <dig>
we obtain the equation system for the functions u and v: ∂∂βuαβ=bvαβ∂∂αvαβ=bvαβ
the elimination of v gives us the following equation for u:  <dig> ∂∂β∂∂βuαβ+abuαβ= <dig>
with u known, one can determine v from the first equation of the system.
if u
α=β
= u
0
and v
α=β
= v
0
are known, we can supply equation  with the boundary condition: uαβ|α=β=u0x∂∂βuαβ|α=β=bv0x
the solution of this equation, already presented in ref  <cit> , is expressed as:  <dig> uαβ=u0α+b∫αβi02abξ−αβ−ξv0ξdξ+ab∫αβi12abξ−αβ−ξβ−ξξ−αu0ξdξ
by reverting to variables x and t and functions g and s and by denoting their values at t =  <dig> as g
0
and s
0
, we obtain: gxt==g0x−v1te−k1t+k2v1+v2e−k1v2+k2v1v1+v2t∫x−v1tx+v2ti02k1k2ξ−x+v1tx+v2t−ξv1+v2s0ξek2−k1x−ξv1+v2dξ+k1k2v1+v2e−k1v2+k2v1v1+v2t∫x−v1tx+v2ti12k1k2ξ−x+v1tx+v2t−ξv1+v2g0ξ⋅⋅ek2−k1x−ξv1+v2x+v2t−ξξ−x+v1tdξ
<dig> s==s0e−k2t+k1v1+v2e−v1+v2t∫x−v1tx+v2ti0v1+v2)ev1+v2g0dξ+k1k2v1+v2e−v1+v2t∫x−v1tx+v2ti1v1+v2)⋅ev1+v2s0dξ
next, we consider g
0
= g
00
im and s
0
= s
00
im = g
00
im as short pulses along the x axis having the width Δ, which is small compared to the interval l, where imx= <dig> if0<x<Δ; <dig> otherwise,
g
<dig> is a linear concentration of g ends at the initial time-point.
numerical calculations based on various values of v
1
; v
2
; k
1
; k
<dig> show that the formulae  indeed represent the solutions of equation .
it must be noted that the pair of functions: gxt=ek2x+v2tv1+v2−k1x−v1tv1+v2sxt=ek2x+v2tv1+v2−k1x−v1tv1+v <dig> also represent a particular solution of equation .
by using this solution at t =  <dig> in , we confirmed that the downstream substitution of the expressions obtained in  results in equality.
since the expression  cannot be easily processed numerically, we approximated the functions g and s by the simpler formulae gf and sf.
these formulae represent an expansion of expression , as follows:  <dig> gfxt==ek2−k1v1+v2x+−k1v2+k2v1v1+v2tg00Δk2v1v1+v2v2i02k1k2−x+v1tx+v2tv1+v2+k1k2v1+v2i12k1k2−x+v1tx+v2tv1+v2x+v2t−x+v1t,if−v2t+Δ<x<v1tg00imx−v1te−k1t,ifv1t≤x≤v1t+Δ <dig> otherwisesfxt==ek2−k1v1+v2x+−k1v2+k2v1v1+v2tg00Δk1v1+v2i02k1k2−x+v1tx+v2tv1+v2+k1k2v1+v2v1v2i12k1k2−x+v1tx+v2tv1+v2−x+v1tx+v2t,if−v2t+Δ<x<v1tg00v1v2imx+v2te−k2t,if−v2t≤x≤−v2t+Δ <dig> otherwise
we are now trying to determine the evolution of short pulses initially situated within a small Δ interval.
a short pulse is concentration of growing  or shrinking  ends of x length in a small Δ interval at a given time.
the analysis for short pulses located close to the origin according to equation  is shown in fig.
<dig>  the short pulses become exponentially weaker with time: the short pulse g moves right at a velocity v
<dig> and its height decreases with time at a rate constant k
<dig>  the short pulse s moves left at a velocity v
<dig>  and its height decreases exponentially with time at a rate constant k
<dig>  we note that the asymptotic approximation we used matches the exact solution fairly well.
it only differs from the exact solution at the top of the short pulse.
however, taking into account that the short pulse width Δ is small, this does not result in a significant distortion of the solution.fig.
<dig> the behavior of the solution of  and the approximate solution of  at t =  <dig> are depicted as dashed and continuous lines, respectively; black and red limes refer to growing and shrinking mts, respectively.
the parameters used are v1 =  <dig>  μm/s, v2 =  <dig>  μm/s, k1 =  <dig>  1/s, k2 =  <dig>  1/s, g00 =  <dig>  the length of the photobleached area is l = 3 μm, and the short pulse width Δ =  <dig>  μm.
arrows denote the direction of the g and s front movements
in conclusion, we found the solution of the system of continuous equations .
an important contribution of this solution is the possibility to differentiate between labeled and unlabeled s ends.
as can be seen in fig.
<dig>  the s ends on the left of the x =  <dig> coordinate point are unlabeled, whereas those on the right are labeled.
solution for finite x coordinate intervals
as mentioned, we found a solution for  for evolution of concentration "pulse" Δ of g and s ends initially located at x =  <dig>  in this case, after time t, g and s ends at the point r located to the right of x =  <dig> gave a marked tubulin track with length r, while the ends on the left did not leave marked tracks.
thus, a positive r can be considered as the length of a marked track.
in this section we are constructing the solution where g and s ends are initially distributed over the 0-l interval of x coordinate  and are then moving in different directions according to .
this can be done if the solution  is convolved with a rectangle over x  =  <dig> if 0 < x < l and rect =  <dig> otherwise).
since for the following analysis we need only the degrading s ends that have no tracks, the auxiliary integration variable xx must be negative and vary from -l to  <dig>  however, when integrating over the initial position of the impulse at 0-Δ, we need to move the limits of xx to -l + Δ, and consider that the Δ interval introduces a small error .
the concentration s
tail
of s ends without tracks  at the coordinate point x and at the time t is: stailxt=1l∫−l+ΔΔrectxsfxx+x,tdxx
the total number of degrading s ends without a tubulin label could be determined by the integration of s
tail
over x. with t =  <dig>  the total number of s ends is g
00
Δv
1
/v
<dig> ; therefore, the fraction s
n
of s ends  among all s ends is: snt=v2v1g00Δ∫0lstailxtdx
this expression permits construction of a solution for a stationary spindle.
frap in a stationary spindle
the essential changes in the mitotic spindle structure occur during prometaphase and at the metaphase-to-anaphase transition.
during metaphase the spindle shape is relatively stable, so that the stationary state appears a good approximation.
in frap experiments on mitotic spindles containing fluorescent tubulin, a rectangular area perpendicular to the spindle axis is photobleached, and fluorescence recovery in this area is recorded  <cit> .
let us consider a stationary case  where: gx=g0e−αxsx=g0v1v2e−αx
in this case, the concentration of g ends in the photobleached area is constant with all g ends incorporating marked tubulin dimers.
at any given time, the concentration of s ends is also constant, but s ends can be of two types: s
n ends that degrade releasing of non-marked photobleached tubulin, and s
m ends that shrink after a short growth phase and therefore degrade releasing fluorescent tubulin.
the fluorescence recovery detection area is usually small compared to the spindle length, and we can neglect concentration changes within this zone.
therefore, we can use g
<dig> and s
0
= g
0
v
1
/v
<dig> instead of g and s; because the spindle is stationary, g
<dig> and s
<dig> do not depend on time.
if large intervals are considered, the equation would need only to account for the exponential dependence of g and s over coordinate x.
growing g
<dig> ends incorporate marked tubulin at a speed of v
<dig>  while the fraction of degrading s
m
ends release marked tubulin at a speed v
<dig>  if m is the length of marked mts within the bleached area, its time derivative would be: ddtmt=v1g0−v2smt
by introducing the fraction of marked tubulin lb = m/e, where e is the total length of marked mts within the zone before photobleaching, and by accounting for g
0
v
1
= s
0
v
<dig> and s
m
= s
0
), we obtain: ddtlbt=v2s0esnt
integration leads to: lbt−lbt0=v2s0e∫t0tsntdt
before photobleaching, most g and s ends within the l interval have marked tubulin at their ends and are related by g
0
v
1
= s
m
v
<dig>  after fluorescence recovery , the g
0
v
1
= s
m
v
<dig> condition is again present.
in such cases, the track length for both g
<dig> and s
m ends before and after fluorescence recovery would be l/ <dig>  fluorescence recovery normalization to the level of fluorescence before bleaching is the accepted mode of experimental data analysis.
in our model, the fluorescence levels before bleaching and after full recovery coincide.
thus, we can normalize the theoretical dependence fluorescence level after full recovery:  <dig> lbt−lbt0lb∞−lbt0=∫t0tsntdt∫t0∞sntdt
is the final dependence that can be used to fit the experimental spindle frap curves.
it is clear that  does not depend on v
2
s
0
/e.
numerical implementation
a pre-anaphase-b frap curve  for drosophila syncytial embryo mitosis was obtained by cheerambathur et al.
<cit> , who also proposed a theoretical model for the frap curve and determined the dynamic parameters of the g and s ends.
approximate parameter values are v
1
= v
2
=  <dig>  μm/s, k
1
=  <dig>  1/s and k
2
=  <dig>  1/s.
the width of the photobleached area near the equator was  <dig>  μm.
they  <cit>  also reported the speed of eb <dig> labeled ends  in their experimental system.
in our model, like in the theoretical model of  <cit> ), v
<dig> and v
<dig> are the speeds in the photobleached area, which is slowly moving to the pole at the speed of the flux .
thus, we calculated  <dig> - <dig>  μm/s for v
<dig>  considered k
1
=  <dig>  1/s, and performed minimization of the root-mean-square deviation of our theoretical curve  from the experimental one  <cit> .
the simple gradient descent method was used to find a local minimum near the dynamic parameters of  <cit> .
the theoretical and experimental curves are shown in fig.
4; the optimal parameters are v
1
=  <dig>  μm/s, v
2
=  <dig>  μm/s, k
1
=  <dig>  1/s, and k
2
=  <dig>  1/s.fig.
<dig> experimental curve for pre-anaphase b  <cit>  and its theoretical approximation .
ex are the corrected experimental points from  <cit> , i.e., the fluorescence level detected just after photobleaching was subtracted from the experimental values; experimental points were then normalized to the fluorescence level before photobleaching.
the theoretical photobleaching curve  was calculated for various values of v <dig>  v <dig>  k <dig>  and k <dig> using mathcad- <dig> software
in our model, we assume that full fluorescence recovery occurs within the bleached area, regardless its position along the spindle.
however, there is experimental evidence that fluorescence recovery is higher near the spindle equator than at the poles  <cit> .
it should be also mentioned that in the stationary case, g and s ends are distributed over the half-spindle in an exponential manner , with the maximum near the pole and the minimum near the equator.
thus, photobleaching near the poles will not only affect mt ends situated in the bleached area but also the ends of “mt fragments” detached from the poles; fluorescence of these “mt fragments” will not be restored during the experiment, lowering the overall level of fluorescence recovery.
if the bleached area is near the equator, this effect would be small.
discussion
the frap technique, although widely used for studying the dynamics of spindle mt behavior, has still some limitations.
these limitations are due to the fact that the kinetic constant value measured in the experiments  has no direct link to the dynamic parameters of mts.
the work of cheerambathur et al.
<cit>  partially addressed this issue by linking this parameter with k
<dig> and k
<dig>  but the authors have not solved the problem completely.
this is because modeling of the fluorescence recovery curve proposed earlier by brust-mascher et al.
<cit>  requires the mechanical equations for the entire spindle to be solved, which in turn requires the knowledge of multiple mechanical parameters.
here, we provided an alternative model for the interpretation of the fluorescence recovery curve.
previously published reports describing mt polymerization and depolymerization dynamics have served as the starting point for our analysis.
based on these studies, we analytically solved the pde, which describes the dynamics and transitions between the states of mt ends.
we then used this solution to solve the equation describing fluorescence recovery in a steady-state spindle and, finally, we obtained a theoretical dependence for the fluorescence recovery curve.
our study provides a model for the frap recovery curve if all four parameters are known.
we also made a model calculation to determine how lowering of each parameter could affect t
1/ <dig> frap time.
as shown in table  <dig>  lowering of v
<dig> or k
<dig> leads to an increase in t
1/ <dig>  while decrease of v
<dig> or k
<dig> decreases t
1/ <dig> table  <dig> alterations of mt growth/shrinking parameters affect t1/ <dig> frap times for half spindles
31
<dig>
33
<dig>
<dig>
<dig>
36
the lower parameter is underlined.
increased and decreased frap times are indicated by bold large numbers and small numbers, respectively
some mutant proteins decrease the rate of the process in which they are involved compared to their wild type counterparts, so we made an attempt to analyze published frap data for three mutant drosophila mitotic proteins.
mini-spindles , a member of the xmap215/tog protein family, concentrates not only at the centrosomes but also at the mt plus-end  <cit> .
based on their own data and pre-existing data buster et al.
<cit>  concluded that msps positively regulates transition from pausing to mt growth state.
although our model does not include mt pausing, a decreased transition to growth means that the constant k
<dig> is decreased.
table  <dig> shows that the decrease in k
<dig> results in an increase in the frap time, which is what has been in fact observed by buster et al.
<cit> .
a second example is concerned with the mast/orbit protein, a clasp family protein that is found at microtubule plus-ends near the kinetochore  <cit> .
rnai-mediated depletion of mast/orbit leads to mt flux inhibition  <cit> , consistent with the finding that mast is involved in the control of microtubule polymerization  <cit> .
thus, mast affects v
<dig> speed in our model and a mutation in mast would decrease v
<dig>  the v
<dig> decrease  decreases frap time, according to the observations of buster et al.
<cit> .
in a third example we consider the eb <dig> protein, a well-known microtubule plus-end marker that increases mt rescue frequency while decreasing pause  <cit> .
loss of eb <dig> would therefore decrease the k
<dig> parameter, which would lead to an increase in frap time as has been in fact observed  <cit> .
the empirical frap curve of salmon et al.
<cit>  permits determination of the dynamic parameter commonly called "tubulin turnover".
here, we have shown that the spindle frap curve depends on four kinetic parameters of mt end polymerization and depolymerization.
using the results of cheerambathur et al.
<cit> , we showed that the solution we found can adequately approximate our experimental frap curve for metaphase spindles.
in addition, we posited that differences in fluorescence recovery between the poles and the equator of metaphase spindles, which are commonly observed in frap experiments, could be explained by the exponential distribution of g and s ends in the half-spindle as proposed in our model.
CONCLUSIONS
here, we provided an alternative model for the interpretation of the fluorescence recovery curve.
based on the previously published reports describing mt polymerization and depolymerization dynamics, we analytically solved the pde, which describes the dynamics and transitions between the states of mt ends.
we then used this solution to solve the equation describing fluorescence recovery in a steady-state spindle and, finally, we obtained a theoretical dependence for the fluorescence recovery curve.
a numerical expression that describes the distribution of polymerizing and depolymerizing microtubule ends as a function of four kinetic parameters is presented.
we also made a model calculation to determine how lowering of each parameter could affect t
1/ <dig> frap time.
as shown in this study, the lowering of v
<dig> or k
<dig> leads to an increase in t
1/ <dig>  while the decrease of v
<dig> or k
<dig> decreases t
1/ <dig>  the numerical expression we elaborated can replace the empirical frap curve described earlier for a spindle comprising fluorescently marked microtubules.
the curve we generated fits well the experimental data.
abbreviations
frapfluorescence recovery after photobleaching
gfpgreen fluorescent protein
mtmicrotubule
pdepartial differential equation

 fluorescence recovery after photobleaching was first introduced in and is a widely used method to study turnover, transport, diffusion and interaction among biological molecules in living specimens. the use of frap has been facilitated by the current availability of microscopes equipped with a laser scanning device. the emergence of fluorescent protein labeling with the green fluorescent protein and its spectral variants has greatly enhanced frap application. in a typical frap experiment, a gfp-labeled structure is rapidly and irreversibly photobleached with a high intensity laser, and fluorescence recovery is recorded as a function of time. the fluorescent molecules then diffuse into the irradiated region, while the non-fluorescent ones diffuse into the unbleached area until equilibrium is reached. the analysis of fluorescence recovery curves yields the diffusion coefficient and the fraction of free, transiently bound and immobilized molecules. for a quantitative description of fluorescence recovery dynamics in frap experiments, several theoretical models have been proposed. these include the pure-diffusion dominant model that considers the recovery rate for weakly bound or free fluorescent molecules, and is defined exclusively by their diffusion; the effective diffusion model, which describes the recovery kinetics of fluorescent molecules that bind tightly the bleached structure, and is also largely defined by diffusion; the reaction dominant model, where diffusion is very fast and molecules rapidly equilibrate after the bleach; and the diffusion phase-binding phase approximation used whenever the contributions of diffusion and binding are coupled. in another early study, salmon et al. used the frap technique to trace the behavior of fluorescently labeled bovine tubulin injected in sea urchin eggs undergoing the first mitotic division. the authors provided evidence that the use of exogenous tubulin accurately represents the in vivo situation, and showed that frap data were best fitted by a negative-going exponential function, and that fluorescence recovery had a half-time of approximately 20 s. taking into account that tubulin dimers diffused back into the bleached area within 1 s, the authors concluded that diffusion was not a limiting factor and that the reaction dominant model was correctly interpreting their frap results. one of the first theoretical description of microtubule growth and degradation was based on data on microtubule dynamics obtained from in vitro experiments. several kinetic models describing microtubule growth were proposed, and one of them, usually referred as the hill’s model, is still widely used. according to this model, mt growth begins with the recruitment of tubulin dimers at a centrosome-associated nucleation site, followed by the polymerization of additional dimers at a rate constant j or depolymerization at a rate constant j. at any length, the mts may switch from the growth mode to a shrinking mode with a constant rate k or from a shrinking mode to a growth mode with a constant rate k a graphical representation of the hill’s model is reported in fig.  this kinetic model was later adapted to describe the polymerization/depolymerization kinetics of microtubules in xenopus egg extracts and to analyze how cyclin a and cyclin b could affect this process. in one of the kinetic regimes, the experimental data were nicely fitting the model, and all of the constants could be defined. in addition to the bounded state where mts are on average disassembling and j 1 k 2 - j 2 k is negative, the authors also described an unbounded state where mts are on average growing.fig. kinetic model of mt growth adapted from hill,. gi and si are the concentrations of growing and shrinking microtubules containing i number of tubulin dimers. c indicates mts with zero length. j and j are the kinetic constants of growth and shrinking, respectively; k k are the constants for rates of growth-to-shrinking and shrinking-to-growth transitions. growing and shrinking mts are depicted as bars, and tubulin dimers are shown as shaded bars the hill’s model was applied to frap-based studies on pre-anaphase b/metaphase spindles of drosophila syncytial embryos expressing gfp-tubulin. the authors concluded that both bounded and unbounded regimes are inadequate to describe the observed frap dynamics. data analysis indicated that tubulin dimers turn over almost entirely during a single cycle of mt shortening and growth, and consequently the recovery time does not depend on the size or position of the bleached region along the metaphase spindle. the recovery halftime after photobleaching is then simply the halftime of this cycle. this critical analysis served as a starting point for us to find other solutions to the kinetic scheme shown in fig.  numerical modeling of a frap experiment shows that the parameters that describe mt turnover can be determined from the modeling data. however, this approach is based on the mechanical model of brust-mascher et al. and requires equations to be solved for the entire spindle, necessitating the knowledge of multiple mechanical parameters. as mentioned above, another approach to analyze mt dynamics is based on chemical kinetics description. this approach operates with a smaller number of dynamic parameters and was verified by in vitro microtubule growth experiments. one of the objectives of the present study is developing a chemical-kinetic description for frap experiments on metaphase spindles. another objective is the analysis of the dependence of spindle frap on parameters describing microtubule end dynamics. frap in a steady-state spindle has a simple relationship with the number of mt growing and shrinking ends within the photobleached area. the dynamics of these ends both with and without specific marks are not easy to evaluate, but can be inferred from the solution of system for elementary intervals. thus, our aim is describing the limit transition from the discrete hill’s equations to a continuous equation, and understanding how the kinetic constants of the hill equation compare to those in the continuous equations. we found a solution for the cauchy problem for partial differential equation for elementary intervals and used this expression to solve the equation for fluorescence recovery in a steady state spindle. finally, we applied our model to the experimental data and performed a quantitative analysis of frap recovery time.  limit passage to the continuous model since we were not able to solve discrete equations describing mt behavior with time, we are using a continuous model. we denote the mt concentration found in the growth phase as g i, where i is the number of polymerized tubulin dimers. similarly, the concentration of microtubules in the shrinking phase is denoted as s i. then, the kinetic equations for the concentration of tubulin dimer appear as: ∂∂tgit=−j1git−gi−1t−k1git+k2sit∂∂tsit=−j2sit−si+1t−k2sit+k1git we can rewrite equation as follows, considering Δx is a small coordinate step: ∂∂tgit=−j1ΔxΔxgit−gi−1t−k1git+k2sit∂∂tsit=−j2ΔxΔxsit−si+1t−k2sit+k1git the values v 1  = j 1 Δx and v 2  = j 2 Δx have a dimension of cm/sec, which correspond to the linear polymerization and depolymerization speeds of a single tubulin dimer. if Δx is negligible compared to the microtubule size, the above equations can be transformed into pdes, where t is a time and x is a coordinate: ∂∂tgxt=−v1∂∂xgxt−k1gxt+k2sxt∂∂tsxt=v2∂∂xsxt−k2sxt+k1gxt as reported below, this transformation will help us to find a mathematical solution for mt behavior. stationary regime and bounded state let us solve the equations for the stationary regime, i.e., when the process does not depend on time, namely when t =  in this case, g and s denote the distribution of growing and shrinking ends over the coordinate x, which corresponds to the current length of mts. v1∂∂xgx=−k1gx+k2sxv2∂∂xsx=−k1gx+k2sx the discriminant of this system, d=k1v1−k2v is either positive or zero. then a characteristic equation: λ2−k2v2−k1v1λ= has two roots, λ1= λ2=k2v1−k1v2v1v2=−α substituting the experimental values of k 1; k 2; v 1; v from reference into the expression for λ results in a positive α. the general solution of the system is gx=c1k2v1+c2k2v1e−αxsx=c1k1v1+c2k2v2e−αx where c and c are arbitrary constants. if at x = 0; g = g 0, and upon x tending to infinity g = g ∞, then gx=g∞+g0−g∞e−αxsx=k1k2g∞+g0−g∞v1v2e−αx using the experimentally observed values of dynamic parameters and arbitrarily assigning g 0  =  and g ∞  =  the system can be visualized as shown in fig.  fig. concentrations of growing and shrinking microtubule ends as a function of coordinate x for a model case if x is measured in μm, the gk 1  = sk ratio holds true for a wide range of large x values, whereas in the case of small x values, which is most pertinent to the experimental situation, an alternative ratio should be applied. a negative-going exponential function for microtubule length observed for the discrete form of equations in the stationary regime was discussed earlier and was successfully used to interpret the mt length distribution in experiments in vitro. based on these experimental data and given the finite number of tubulin molecules in the cell, g ∞  >  is not a viable option and so g ∞ should be set as in this case, a simple ratio gv 1  = sv is true and corresponds to a situation where the spindle remains unchanged with time, and the numbers of growing and shrinking microtubules are equal. analytical solution of the equations for elementary intervals the system of continuous equations was first published by dogterom and leibler; the authors stated that it could be solved analytically with the border conditions set at t =  and with all mts having zero length. however, they did not publish the solution. in the following section we find this solution. the solution refers to a frap experiment that involves photobleaching of fluorescently labeled mts in a rectangular region of width l perpendicular to the spindle axis, followed by the analysis of the dynamics of fluorescence recovery in the bleached region. in a real experiment both g and s mt ends are distributed uniformly across this region at the initial time-point after bleaching. let us consider the elementary interval Δ of l, which initially contains g and s unmarked ends. the g ends will start to produce a marked tubulin track and will move towards the spindle equator from the elementary interval, while the s ends will move to the opposite direction. some g ends will then turn into s ends and some s ends will turn into g ends. to describe the behavior of g and s ends it is necessary to solve the cauchy problem for. we thus substituted the x and t variables with α and β, and converted the functions g and s into u and v: α=x−v1tβ=x+v2tgxt=uαβe−αβ+basxt=vαβe−αβ+baa=k1v1+v2b=k2v1+v we obtain the equation system for the functions u and v: ∂∂βuαβ=bvαβ∂∂αvαβ=bvαβ the elimination of v gives us the following equation for u: ∂∂β∂∂βuαβ+abuαβ= with u known, one can determine v from the first equation of the system. if u α=β  = u 0 and v α=β  = v 0 are known, we can supply equation with the boundary condition: uαβ|α=β=u0x∂∂βuαβ|α=β=bv0x the solution of this equation, already presented in ref, is expressed as: uαβ=u0α+b∫αβi02abξ−αβ−ξv0ξdξ+ab∫αβi12abξ−αβ−ξβ−ξξ−αu0ξdξ by reverting to variables x and t and functions g and s and by denoting their values at t =  as g 0 and s 0, we obtain: gxt==g0x−v1te−k1t+k2v1+v2e−k1v2+k2v1v1+v2t∫x−v1tx+v2ti02k1k2ξ−x+v1tx+v2t−ξv1+v2s0ξek2−k1x−ξv1+v2dξ+k1k2v1+v2e−k1v2+k2v1v1+v2t∫x−v1tx+v2ti12k1k2ξ−x+v1tx+v2t−ξv1+v2g0ξ⋅⋅ek2−k1x−ξv1+v2x+v2t−ξξ−x+v1tdξ s==s0e−k2t+k1v1+v2e−v1+v2t∫x−v1tx+v2ti0v1+v2)ev1+v2g0dξ+k1k2v1+v2e−v1+v2t∫x−v1tx+v2ti1v1+v2)⋅ev1+v2s0dξ next, we consider g 0  = g 00 im and s 0  = s 00 im = g 00 im as short pulses along the x axis having the width Δ, which is small compared to the interval l, where imx= if0<x<Δ; otherwise, g is a linear concentration of g ends at the initial time-point. numerical calculations based on various values of v 1; v 2; k 1; k show that the formulae indeed represent the solutions of equation. it must be noted that the pair of functions: gxt=ek2x+v2tv1+v2−k1x−v1tv1+v2sxt=ek2x+v2tv1+v2−k1x−v1tv1+v also represent a particular solution of equation. by using this solution at t =  in, we confirmed that the downstream substitution of the expressions obtained in results in equality. since the expression cannot be easily processed numerically, we approximated the functions g and s by the simpler formulae gf and sf. these formulae represent an expansion of expression, as follows: gfxt==ek2−k1v1+v2x+−k1v2+k2v1v1+v2tg00Δk2v1v1+v2v2i02k1k2−x+v1tx+v2tv1+v2+k1k2v1+v2i12k1k2−x+v1tx+v2tv1+v2x+v2t−x+v1t,if−v2t+Δ<x<v1tg00imx−v1te−k1t,ifv1t≤x≤v1t+Δ otherwisesfxt==ek2−k1v1+v2x+−k1v2+k2v1v1+v2tg00Δk1v1+v2i02k1k2−x+v1tx+v2tv1+v2+k1k2v1+v2v1v2i12k1k2−x+v1tx+v2tv1+v2−x+v1tx+v2t,if−v2t+Δ<x<v1tg00v1v2imx+v2te−k2t,if−v2t≤x≤−v2t+Δ otherwise we are now trying to determine the evolution of short pulses initially situated within a small Δ interval. a short pulse is concentration of growing or shrinking ends of x length in a small Δ interval at a given time. the analysis for short pulses located close to the origin according to equation is shown in fig.  the short pulses become exponentially weaker with time: the short pulse g moves right at a velocity v and its height decreases with time at a rate constant k the short pulse s moves left at a velocity v and its height decreases exponentially with time at a rate constant k we note that the asymptotic approximation we used matches the exact solution fairly well. it only differs from the exact solution at the top of the short pulse. however, taking into account that the short pulse width Δ is small, this does not result in a significant distortion of the solution.fig. the behavior of the solution of and the approximate solution of at t =  are depicted as dashed and continuous lines, respectively; black and red limes refer to growing and shrinking mts, respectively. the parameters used are v1 =   μm/s, v2 =   μm/s, k1 =  1/s, k2 =  1/s, g00 =  the length of the photobleached area is l = 3 μm, and the short pulse width Δ =   μm. arrows denote the direction of the g and s front movements in conclusion, we found the solution of the system of continuous equations. an important contribution of this solution is the possibility to differentiate between labeled and unlabeled s ends. as can be seen in fig.  the s ends on the left of the x =  coordinate point are unlabeled, whereas those on the right are labeled. solution for finite x coordinate intervals as mentioned, we found a solution for for evolution of concentration "pulse" Δ of g and s ends initially located at x =  in this case, after time t, g and s ends at the point r located to the right of x =  gave a marked tubulin track with length r, while the ends on the left did not leave marked tracks. thus, a positive r can be considered as the length of a marked track. in this section we are constructing the solution where g and s ends are initially distributed over the 0-l interval of x coordinate and are then moving in different directions according to. this can be done if the solution is convolved with a rectangle over x  =  if 0 < x < l and rect =  otherwise). since for the following analysis we need only the degrading s ends that have no tracks, the auxiliary integration variable xx must be negative and vary from -l to however, when integrating over the initial position of the impulse at 0-Δ, we need to move the limits of xx to -l + Δ, and consider that the Δ interval introduces a small error. the concentration s tail of s ends without tracks at the coordinate point x and at the time t is: stailxt=1l∫−l+ΔΔrectxsfxx+x,tdxx the total number of degrading s ends without a tubulin label could be determined by the integration of s tail over x. with t =  the total number of s ends is g 00 Δv 1 /v; therefore, the fraction s n of s ends among all s ends is: snt=v2v1g00Δ∫0lstailxtdx this expression permits construction of a solution for a stationary spindle. frap in a stationary spindle the essential changes in the mitotic spindle structure occur during prometaphase and at the metaphase-to-anaphase transition. during metaphase the spindle shape is relatively stable, so that the stationary state appears a good approximation. in frap experiments on mitotic spindles containing fluorescent tubulin, a rectangular area perpendicular to the spindle axis is photobleached, and fluorescence recovery in this area is recorded. let us consider a stationary case where: gx=g0e−αxsx=g0v1v2e−αx in this case, the concentration of g ends in the photobleached area is constant with all g ends incorporating marked tubulin dimers. at any given time, the concentration of s ends is also constant, but s ends can be of two types: s n ends that degrade releasing of non-marked photobleached tubulin, and s m ends that shrink after a short growth phase and therefore degrade releasing fluorescent tubulin. the fluorescence recovery detection area is usually small compared to the spindle length, and we can neglect concentration changes within this zone. therefore, we can use g and s 0  = g 0 v 1 /v instead of g and s; because the spindle is stationary, g and s do not depend on time. if large intervals are considered, the equation would need only to account for the exponential dependence of g and s over coordinate x. growing g ends incorporate marked tubulin at a speed of v while the fraction of degrading s m ends release marked tubulin at a speed v if m is the length of marked mts within the bleached area, its time derivative would be: ddtmt=v1g0−v2smt by introducing the fraction of marked tubulin lb = m/e, where e is the total length of marked mts within the zone before photobleaching, and by accounting for g 0 v 1  = s 0 v and s m  = s 0 ), we obtain: ddtlbt=v2s0esnt integration leads to: lbt−lbt0=v2s0e∫t0tsntdt before photobleaching, most g and s ends within the l interval have marked tubulin at their ends and are related by g 0 v 1  = s m v after fluorescence recovery, the g 0 v 1  = s m v condition is again present. in such cases, the track length for both g and s m ends before and after fluorescence recovery would be l/ fluorescence recovery normalization to the level of fluorescence before bleaching is the accepted mode of experimental data analysis. in our model, the fluorescence levels before bleaching and after full recovery coincide. thus, we can normalize the theoretical dependence fluorescence level after full recovery: lbt−lbt0lb∞−lbt0=∫t0tsntdt∫t0∞sntdt is the final dependence that can be used to fit the experimental spindle frap curves. it is clear that does not depend on v 2 s 0 /e. numerical implementation a pre-anaphase-b frap curve for drosophila syncytial embryo mitosis was obtained by cheerambathur et al., who also proposed a theoretical model for the frap curve and determined the dynamic parameters of the g and s ends. approximate parameter values are v 1  = v 2  =   μm/s, k 1  =  1/s and k 2  =  1/s. the width of the photobleached area near the equator was  μm. they also reported the speed of eb labeled ends in their experimental system. in our model, like in the theoretical model of ), v and v are the speeds in the photobleached area, which is slowly moving to the pole at the speed of the flux. thus, we calculated -  μm/s for v considered k 1  =  1/s, and performed minimization of the root-mean-square deviation of our theoretical curve from the experimental one. the simple gradient descent method was used to find a local minimum near the dynamic parameters of. the theoretical and experimental curves are shown in fig. 4; the optimal parameters are v 1  =   μm/s, v 2  =   μm/s, k 1  =  1/s, and k 2  =  1/s.fig. experimental curve for pre-anaphase b and its theoretical approximation. ex are the corrected experimental points from, i.e., the fluorescence level detected just after photobleaching was subtracted from the experimental values; experimental points were then normalized to the fluorescence level before photobleaching. the theoretical photobleaching curve was calculated for various values of v v k and k using mathcad- software in our model, we assume that full fluorescence recovery occurs within the bleached area, regardless its position along the spindle. however, there is experimental evidence that fluorescence recovery is higher near the spindle equator than at the poles. it should be also mentioned that in the stationary case, g and s ends are distributed over the half-spindle in an exponential manner, with the maximum near the pole and the minimum near the equator. thus, photobleaching near the poles will not only affect mt ends situated in the bleached area but also the ends of “mt fragments” detached from the poles; fluorescence of these “mt fragments” will not be restored during the experiment, lowering the overall level of fluorescence recovery. if the bleached area is near the equator, this effect would be small. discussion the frap technique, although widely used for studying the dynamics of spindle mt behavior, has still some limitations. these limitations are due to the fact that the kinetic constant value measured in the experiments has no direct link to the dynamic parameters of mts. the work of cheerambathur et al. partially addressed this issue by linking this parameter with k and k but the authors have not solved the problem completely. this is because modeling of the fluorescence recovery curve proposed earlier by brust-mascher et al. requires the mechanical equations for the entire spindle to be solved, which in turn requires the knowledge of multiple mechanical parameters. here, we provided an alternative model for the interpretation of the fluorescence recovery curve. previously published reports describing mt polymerization and depolymerization dynamics have served as the starting point for our analysis. based on these studies, we analytically solved the pde, which describes the dynamics and transitions between the states of mt ends. we then used this solution to solve the equation describing fluorescence recovery in a steady-state spindle and, finally, we obtained a theoretical dependence for the fluorescence recovery curve. our study provides a model for the frap recovery curve if all four parameters are known. we also made a model calculation to determine how lowering of each parameter could affect t 1/ frap time. as shown in table  lowering of v or k leads to an increase in t 1/ while decrease of v or k decreases t 1/ table alterations of mt growth/shrinking parameters affect t1/ frap times for half spindles 31 33 36 the lower parameter is underlined. increased and decreased frap times are indicated by bold large numbers and small numbers, respectively some mutant proteins decrease the rate of the process in which they are involved compared to their wild type counterparts, so we made an attempt to analyze published frap data for three mutant drosophila mitotic proteins. mini-spindles, a member of the xmap215/tog protein family, concentrates not only at the centrosomes but also at the mt plus-end. based on their own data and pre-existing data buster et al. concluded that msps positively regulates transition from pausing to mt growth state. although our model does not include mt pausing, a decreased transition to growth means that the constant k is decreased. table  shows that the decrease in k results in an increase in the frap time, which is what has been in fact observed by buster et al.. a second example is concerned with the mast/orbit protein, a clasp family protein that is found at microtubule plus-ends near the kinetochore. rnai-mediated depletion of mast/orbit leads to mt flux inhibition, consistent with the finding that mast is involved in the control of microtubule polymerization. thus, mast affects v speed in our model and a mutation in mast would decrease v the v decrease decreases frap time, according to the observations of buster et al.. in a third example we consider the eb protein, a well-known microtubule plus-end marker that increases mt rescue frequency while decreasing pause. loss of eb would therefore decrease the k parameter, which would lead to an increase in frap time as has been in fact observed. the empirical frap curve of salmon et al. permits determination of the dynamic parameter commonly called "tubulin turnover". here, we have shown that the spindle frap curve depends on four kinetic parameters of mt end polymerization and depolymerization. using the results of cheerambathur et al., we showed that the solution we found can adequately approximate our experimental frap curve for metaphase spindles. in addition, we posited that differences in fluorescence recovery between the poles and the equator of metaphase spindles, which are commonly observed in frap experiments, could be explained by the exponential distribution of g and s ends in the half-spindle as proposed in our model.  here, we provided an alternative model for the interpretation of the fluorescence recovery curve. based on the previously published reports describing mt polymerization and depolymerization dynamics, we analytically solved the pde, which describes the dynamics and transitions between the states of mt ends. we then used this solution to solve the equation describing fluorescence recovery in a steady-state spindle and, finally, we obtained a theoretical dependence for the fluorescence recovery curve. a numerical expression that describes the distribution of polymerizing and depolymerizing microtubule ends as a function of four kinetic parameters is presented. we also made a model calculation to determine how lowering of each parameter could affect t 1/ frap time. as shown in this study, the lowering of v or k leads to an increase in t 1/ while the decrease of v or k decreases t 1/ the numerical expression we elaborated can replace the empirical frap curve described earlier for a spindle comprising fluorescently marked microtubules. the curve we generated fits well the experimental data. abbreviations frapfluorescence recovery after photobleaching gfpgreen fluorescent protein mtmicrotubule pdepartial differential equation 
the empirical frap curve proposed earlier permits determination of only one such dynamic parameter, commonly called the "tubulin turnover". a numerical expression that describes the distribution of polymerizing and depolymerizing microtubule ends as a function of four kinetic parameters is presented. mitotic spindlefrap microtubulesfluorescently marked tubulingrowing/shrinking microtubule endspde the international conference on bioinformatics of genome regulation and structure\systems biology novosibirsk, russia august- september issue-copyright-statement© the author 2017 spindle frap curves depend on the kinetic parameters of microtubule polymerization and depolymerization. the aim of our study was to build a frap curve based on an already known kinetic model of microtubule growth. in addition, a theoretical frap curve for the metaphase spindle is constructed using previously published dynamic parameters. the numerical expression we elaborated can replace the empirical frap curve described earlier for a spindle comprising fluorescently marked microtubules. the curve we generated fits well the experimental data.
 spindle frap curves depend on the kinetic parameters of microtubule polymerization and depolymerization. the empirical frap curve proposed earlier permits determination of only one such dynamic parameter, commonly called the "tubulin turnover". the aim of our study was to build a frap curve based on an already known kinetic model of microtubule growth.  a numerical expression that describes the distribution of polymerizing and depolymerizing microtubule ends as a function of four kinetic parameters is presented. in addition, a theoretical frap curve for the metaphase spindle is constructed using previously published dynamic parameters.  the numerical expression we elaborated can replace the empirical frap curve described earlier for a spindle comprising fluorescently marked microtubules. the curve we generated fits well the experimental data. keywords mitotic spindlefrap microtubulesfluorescently marked tubulingrowing/shrinking microtubule endspde the international conference on bioinformatics of genome regulation and structure\systems biology novosibirsk, russia august- september issue-copyright-statement© the author 2017 