Using Fluorescence Recovery After Photobleaching data to uncover filament dynamics

Fluorescence Recovery After Photobleaching (FRAP) has been extensively used to understand molecular dynamics in cells. This technique when applied to soluble, globular molecules driven by diffusion is easily interpreted and well understood. However, the classical methods of analysis cannot be applied to anisotropic structures subjected to directed transport, such as cytoskeletal filaments or elongated organelles transported along microtubule tracks. A new mathematical approach is needed to analyze FRAP data in this context and determine what information can be obtain from such experiments. To address these questions, we analyze fluorescence intensity profile curves after photobleaching of fluorescently labelled intermediate filaments anterogradely transported along microtubules. We apply the analysis to intermediate filament data to determine information about the filament motion. Our analysis consists of deriving equations for fluorescence intensity profiles and developing a mathematical model for the motion of filaments and simulating the model. Two closed forms for profile curves were derived, one for filaments of constant length and one for filaments with constant velocity, and three types of simulation were carried out. In the first type of simulation, the filaments have random velocities which are constant for the duration of the simulation. In the second type, filaments have random velocities which instantaneously change at random times. In the third type, filaments have random velocities and exhibit pausing between velocity changes. Our analysis shows: the most important distribution governing the shape of the intensity profile curves obtained from filaments is the distribution of the filament velocity. Furthermore, filament length which is constant during the experiment, had little impact on intensity profile curves. Finally, gamma distributions for the filament velocity with pauses give the best fit to asymmetric fluorescence intensity profiles of intermediate filaments observed in FRAP experiments performed in polarized migrating astrocytes. Our analysis also shows that the majority of filaments are stationary. Overall, our data give new insight into the regulation of intermediate filament dynamics during cell migration.


Introduction
Living organisms are in constant dynamic equilibrium. In cells, many structures appear generally static but are, in fact, formed of molecules continuously moving and exchanging with the surrounding. Fluorescence Recovery After Photobleaching (FRAP), developed in the 1970s, is an essential tool for understanding molecular dynamics within a cell [1][2][3]. The typical setup for a FRAP experiment involves a fluorescent probe, a microscope, and some method of photobleaching [4]. A portion of the domain where the molecule of interest is present is bleached and the recovery of fluorescence in that region is imaged over time. In order to gain quantitative information on molecular dynamics, mathematical models of diffusion are typically used. These include models of diffusion in inhomogeneous media [5], models of diffusion and binding using reaction-diffusion equations [6][7][8], and advection-reaction-diffusion models of active transport and diffusion [9]. In some instances, when diffusion parameters are not of interest, simpler ordinary differential equation models are used to elicit information [10]. All these models deal with analysis of soluble, generally globular, molecules. Until now there has been no analysis of FRAP data regarding the dynamics of filamentous structures.
The main example we have in mind is that of short term transport of mature intermediate filaments (IFs), one of three major fibrous structural components of the cytoskeleton. They form a filamentous network spreading throughout the cell cytoplasm and this network together with actin filaments and microtubules, plays a key role in cell polarity and migration [11]. In migrating astrocytes (the type of glial cell used in our migration experiments), the dynamics of the IF network is mainly driven by microtubule and actin mediated transport [12,13]. Deterministic and stochastic mathematical models have been developed to describe the motion of IF driven by antagonistic molecular motors along microtubules [14,15]. In [12], FRAP experiments of IFs were used to better understand how the IF network global dynamics are regulated in migrating and non-migrating glial cells. They showed that, during cell polarization, IF transport is mainly anterograde, oriented from the cell center to the cell periphery, and this bias was due to the inhibition of the retrograde transport of IFs by CDC42-driven polarity signaling. However, due to the high density of the IF network, it was not possible to quantify the dynamics of IFs at the single filament level. Hence, there is a need for a mathematical model to better understand collective IF transport using FRAP data.
There are two types of data gathered when conducting FRAP experiments. The first is a series of time measurements, called the profile curves, which show the profile of fluorescence intensity plotted along the direction of migration (X 1 -direction) across the bleached region, and integrated along the perpendicular direction (X 2 -direction) as depicted in Fig 1a. The second type of data is the total fluorescent intensity of the bleached region as a function of time (fluorescence intensity vs time after bleaching) and is called a fluorescence recovery curve.
Since the former give more information we focus on those.
Our goal here is to determine what information can be obtained from the profile curves obtained from FRAP data. In particular, is it possible to infer the mode of filament motion (constant or variable velocity, with or without pausing), or information about filament velocity and length?

Methods
Experiments were performed in astrocytes [12] which are the major glial cells of the central nervous system. Astrocyte polarization and migration was induced by a scratch wound. We studied the IF dynamics in cells at the wound 1-2 hours after wounding, when cells are polarizing and migrate perpendicularly to the wound axis (see Fig 1a). IF are transported by molecular motors along the polarized microtubule network, which is aligned with the front-to-rear polarity axis, with the microtubule minus ends concentrated near the cell center and the plus-ends growing towards the cell's leading edge [16]. Therefore, IF motion is also mostly parallel to the direction of migration. The bleached region is rectangular and its height is perpendicular to the protrusion, i.e. parallel to the wound (see Fig 1a). In these conditions, the fluorescence profile curve is asymmetric, due to the inhibition of the dynein-dependent retrograde transport of IFs and reflects the anterograde transport of IF dominated by kinesin motors along the polarized network of microtubule [12]. Hence in these experiments most of the transported IFs move from the center to the front of the cell and across the width of the bleached region. This allows us to reduce the problem to one dimension in the direction of the width of the bleached region (X 1 direction in Fig 1a). Furthermore, we assume the density of the filaments is uniform in the direction of the height of the bleached region (X 2 direction in Fig 1a).
The time scale of the experiments is less than 30 minutes and the fluorescence only comes back from the edges of the bleached region. Hence we assume diffusion and remodelling of the filaments due to polymerization/depolymerization, subunit exchange (which occurs on a timescale of hours), fusion or severing are negligible [12,17]. Thus, the length of the filaments is assumed to be fixed during the observation time and the active transport of filaments moving from the cell center to the cell front is the major mechanism causing fluorescence recovery. In this model crowding effects or interactions with other filaments or organelles are not taken into consideration.
We will use mathematical modeling to determine what information can be obtained from FRAP data in the context of directional transport of elongated structures. In particular, we will focus on what characteristics of the velocity, length, and pause of filaments can be deduced from the experimental data. Due to the one directional transport, velocity and speed are synonymous in this work.

The mathematical model
filaments to move in the positive X 1 direction (direction of the width of the bleached region) independent of the other filaments. Formally the mathematical model of motion can be stated as: let x(t) be the position of the right endpoint of the fiber at time t, with v � 0 the velocity, x 0 the initial position, and ℓ the fiber length which does not change in time. Then Fig 1. A schematic of the model and an example of a FRAP experiment. Panel (a) shows the schematic representation of a cell migrating to the right with several cytoplasmic IFs. Selected velocities are denoted with purple arrows under the filaments. The right endpoints of the IFs are denoted x 0 and one is labeled with its length ℓ and velocity v. The right endpoints are distributed between 0 and F. Panels (b)-(g) show data from a FRAP experiment. Panels (b)-(d) show the profile curves for (e)-(g) respectively. The curves in (b)-(d) come from data in a subregion of (e)-(g); for example the yellow box in (f). Panels (e)-(g) show fluorescence images taken from a FRAP experiment performed on vimentin-EGFP expressing astrocytes, located at the edge of a wound, 1h after wounding of the monolayer. The rear-to-front polarity axis similar to the direction of migration is indicated by an arrow in (a), (e) and (h). Data are before (b) and (e), just after (c) and (f), and 2 minutes after bleaching (d) and (g). Panels (h)-(j) are blowups of the region indicated by the yellow box in (f) and show the domain and setup for the simulations. For the clarity, the vertical coordinate of each filament stays the same. The photobleached region (shown in white in (i) and (j)) is between y 0 and y 1 . Several filaments are shown with their initial right endpoint depicted by a dot in (h)-(j). The red filaments are not relevant for our mathematical analysis and simulations since we consider only right moving filaments. Of the remaining filaments, the green ones are unbleached and the blue ones are partially bleached. In (i) the squares denotedx 0 show where the new right endpoints will be after bleaching. In the mathematical analysis and simulations it is as if all the filaments to the right of y 0 are bleached since we only consider right moving filaments.
https://doi.org/10.1371/journal.pcbi.1010573.g001 and the left endpoint is at x(t) − ℓ. We consider the random variables X 0 , V, and L, which we assume to be independent, corresponding to x 0 , v, and ℓ. We define the corresponding random process X(t) = Vt + X 0 (see Fig 1a for a depiction of the setup).
We compare our results with experiments from polarized migrating astrocytes showing asymmetric fluorescence intensity profile curves with reduced retrograde transport [12]. Thus we consider only filaments which move to the right (having non-negative velocities). Following the observations that neurofilaments, a type of IFs observed in neurons, display a stopand-go motion [13,18], we assume that the filaments can have a stop-and-go behavior, moving for a period of time T m at a velocity V and then stopping for a period of time T s . The cycle repeats, with the motion time, velocity, and stop time for each cycle being independent of those for the other cycles.
We now consider the experimental process of FRAPing in the context of IFs (or any elongated objects transported unidirectionally). We assume the right endpoints of the filaments are uniformly distributed on the interval [0, F]. Suppose that y 0 < y 1 (the start and end of the bleach zone) and a region ½y 0 ; y 1 � � X 2 is bleached. We also assume y 0 � F to ensure that there are filaments in the bleached region at the time of bleaching. The bleaching process does not change the underlying behavior of the fibers, but the bleached portions of the fibers are no longer visible and do not contribute to the fluorescence intensity profiles (see Fig 1).
We categorize filaments into three types: unbleached filaments, partially bleached filaments, and entirely bleached filaments. Since we only use fluorescence intensity data, we only consider the unbleached portions of filaments as these are the parts of the filaments that fluoresce. For our analysis we restrict the filaments to have non-negative velocities. We need not consider filaments which do not extend to the left of y 0 (i.e., for which x 0 − ℓ > y 0 ; red filaments in Fig 1h) and filaments whose left endpoint is to right of y 1 (none of these filaments are shown in Fig 1h-1j). The rest of the partially bleached filaments (the blue filaments in Fig 1h) are defined to have a new right endpoint where the unbleached portion of the filament to the left of the bleached portion starts (marked by squares on the filaments in Fig 1i). And of course, the unbleached filaments which lie to the left of the bleached region (the green filaments in Fig  1h-1j) are considered. Mathematically we say the right endpoints of the filaments under consideration are defined in the following manner: x 0 otherwise: (

Simulations
There are five random variables with their associated distributions in the model. The initial setup of the filaments is determined by two random variables: the initial position of the right endpoint X 0 and the fiber length L. The movement of filaments is governed by V the filament velocity, T s the pausing time of the filament, and T m the time the filament is moving. We denote the distributions for all the random variables as follows: • m X 0 governs the initial right endpoint X 0 of the filaments, • μ L governs the length L of each filament, • μ V governs the velocity V of the filament during each period of motion, • μ on governs the duration of time each filament moves before pausing or changing velocity, and • μ off governs the duration of time each filament remains stationary (or pausing) before moving again.
The distribution m X 0 that governs the initial position of the right endpoint of filaments is always a uniform distribution on the interval [0, F], for some F > 0 and y 0 � F. The distributions of the moving and pausing times, μ on and μ off , are also always uniform and only the parameters are unknown. The distribution types and their characteristic parameters for filament length and velocity, μ L and μ V , are also unknown. The types of distributions considered are Dirac delta (or deterministic), uniform, normal, and gamma distributions. Dirac delta distribution would reflect a case where the velocity of filament is precisely set. Uniform distribution of velocities would mean that all velocities are equiprobable, and that there is no internal control of filament velocity within the cell. For example, the velocity could be between 0 and the maximal free load velocity of kinesin motors. The Gaussian distribution would reflect the case where velocities are distributed symmetrically around an average value, suggesting that the control of filament velocity is noisy but symmetrical. The gamma distribution would reflect the fact the velocity is asymmetrically distributed with a higher contribution of slow filaments. The gamma distribution would empirically describe the asymmetric velocity distributions predicted in the transport of cargoes when friction plays an important role [19,20], as it is the case for IF [12]. The Dirac distribution depends on one free parameter and the other three have two free parameters. The mean denoted μ and standard deviation denoted σ are defined in the standard manner for the uniform and Gaussian distributions and for the gamma distribution μ = kθ and s ¼ ffi ffi ffi k p y where k and θ are the free shape and scale parameters respectively. Hence the determination of the appropriate distribution and their relevant parameters to use is the primary objective of this study.
We numerically simulated the FRAP experiments by moving the filaments and calculating how many filaments are in the bleached region. We did this in three different ways: 1) each filament gets a different velocity determined by μ V but it remains constant throughout the simulation; 2) after T on time units have elapsed, where T on is determined by μ on , velocities change but always come from the same distribution, μ V ; and 3) filaments have velocities which change, again determined by μ V , but they stop in between velocity changes. Thus, there is a cycle for each filament of duration T on + T off where T on is a μ on distributed random variable and is the time the filament is moving (the motion is on) and T off is a μ off distributed random variable and is the time the filament is stationary (the motion is off). We refer to these simulations as type 1, 2, and 3 simulations. Fig 2 depicts three filaments for each type of simulation. Depending on the simulation the length of the filaments is either fixed or uniformly, normally, or gamma distributed. Similarly, the velocity is either fixed or uniformly, normally, or gamma distributed. For the velocity the normal and gamma distributions are truncated so no velocities are greater than 40 microns per minute (and for the normal distribution the velocities are all positive). For type 2 simulations the lengths of time during which velocity is fixed are uniformly distributed with a specified mean τ on . For type 3 simulations the stop and run times are uniformly distributed with specified means τ off and τ on .
We use MATLAB to perform the simulations. This work is driven by experimental data, however the variables of interest (velocity, length, pausing, and moving time) are not observable. Hence we model the variables of interest as random variables with underlying distributions. As previously motivated we use 4 types of distributions. We estimate only the parameters (mean and variance) of these distributions. Hence we fit the simulated fluorescence intensities to FRAP data by calibrating these distributions using the MATLAB function fmincon, a nonlinear optimizer which finds the minimum of a constrained nonlinear multivariable function. The result of our fitting specifies which type of distribution to use and its relevant free parameters for each variable of interest (velocity, pausing, and moving times). We found that information about length is not encoded in the profile curves, thus, to reduce the model complexity during the fitting process, we fixed the length to be a uniformly distributed random variable with μ = 5.025.

Cell culture.
Primary rat astrocytes were prepared as previously described [21] according to the guidelines approved by the French Ministry of Agriculture and following European standards. For scratch-induced migration assays, cells were seeded on poly-L-ornithine-precoated coverslips for immunofluorescence or 35-mm glass-bottomed culture dishes (MatTek Corporation) for videomicroscopy. Cells were grown to confluence in DMEM with 1 g/l glucose and supplemented with 10% FBS (Invitrogen), 1% penicillin-streptomycin (Thermo Fisher Scientific), and 1% amphotericin B (Thermo Fisher Scientific). On the day of the experiment, cells were scratched with a blunt-ended microinjection needle, creating a 300μm-wide wound to trigger cell migration.

Cell transfection.
Starting from a 10 cm diameter petri dish, primary astrocytes grown to confluence were trypsinized and electroporated with a Nucleofector machine (Lonza)  (4) where all filaments have the same fixed velocity but random (time-independent) lengths. Panel (e) the particular situation/case of type 1 simulations theoretically resolved in Eq (5) where all the filaments have the same length but different (time-independent) velocities that are randomly selected.
https://doi.org/10.1371/journal.pcbi.1010573.g002 using 5 μg of vimentin-EGFP DNA. We have previously shown that EGFP-tagged vimentin co-polymerizes with the endogenous IF proteins and fluorescently labels the whole astrocytic IF network. Therefore labeling vimentin fluorescently is enough to follow the dynamics of the complete/whole IF network [12]. Medium was changed the day after transfection.

Live-cell imaging.
Nucleofected primary astrocytes were seeded on 35-mm glassbottomed dishes and grown to confluence for 4 days. On the day before wounding, the medium was changed to a phenol red-free DMEM supplemented with 10% serum. The monolayer was wounded and cells were monitored between 1 and 2 hours after wounding, allowing them to grow a polarized protrusion [22]. Videos were acquired on a spinning-disk confocal microscope (PerkinElmer) equipped with an electron-multiplying charge-coupled device camera and either a 63×, 1.4 NA objective or a 100×, 1.4 NA objective.

Results
In this section we give the results of the mathematical theory, the three types of simulations, and the experimental data. We divide it into five main subsections: first, we explain the theoretical results; second, we consider what can be learned from the initial setup; third, we compare type 1 simulations (where the velocity for each filament is fixed for the duration of the simulation but each filament's velocity can be different) with the theoretical results derived in Eqs (4) and (5); fourth, we compare results from type 1 simulations (where each filament can have a different but fixed velocity) with type 2 simulations (where the velocity can abruptly change to a new value during the simulation) and with type 3 simulations (where the filament pauses before changing velocity); and finally, we compare the theory and results from simulations of type 1 and type 3 with experimental data.

Theoretical results
Based on the filament motion model assumed in this work and the description of the experimental setup described above, we are now deriving closed forms for the profile curves. Two simplifications allow the derivation of two equations for the profile curves valid under the corresponding assumptions. First, we assume that all the filaments have the same fixed velocity (a special case of type 1 simulations where all the filaments have a fixed velocity which is the same, i.e., a Dirac delta distribution which gives all the filaments the same velocity, see Fig 2d). Thus the velocity is deterministic and no longer random. We then derive the corresponding profile curves in Eq (4). Second, we used a random non-fixed velocity and we fix an identical filament length for all filaments, see Fig 2e. This allows us to derive Eq (5).
Let y represent an arbitrary point in the bleached zone, i.e., y 2 [y 0 , y 1 ]. The probability that some part of the filament is at y is given by PðfXðtÞ > yg \ fXðtÞ À L < ygÞ: Recall we are interested only in filaments which enter the bleached region from the left. There are two types of filaments, ones which are not bleached and ones which are partially bleached. Thus we define two sets, given y and t. Let be the set of values corresponding to filaments that are not bleached, because their right endpoints are located before the bleached region (the green filaments in Fig 1h-1j). Likewise let Bðt; yÞ ¼ fðx 0 ; v; 'Þjy 0 þ vt > y and y 0 þ vt À ð' À ðx 0 À y 0 ÞÞ < y and x 0 > y 0 and x 0 À ' < y 0 g ð2Þ be the set of values corresponding to the non-bleached region of filaments that are only partially bleached (the blue filaments in Fig 1i-1j). We note that these two sets are disjoint so If we let P be the distribution of (X 0 , L, V), and let E be the corresponding expectation, then the first probability becomes and we can obtain the second probability similarly. The profile curves are scaled versions of while the fluorescence recovery curve is a scaled version of GðtÞ ¼ R y 1 y 0 Hðt; yÞ dy. Recall that it is assumed that filaments right endpoints are uniformly distributed; so let X 0 be uniformly distributed on the interval [0, F] with F � y 0 . Moreover, define In order to compute H(t, y) explicitly, we make one of two simplifying assumptions. First, suppose V = v with probability one (see Fig 2d). In this case it is convenient to work in traveling wave coordinates so we let w ≔ y − vt. In this scenario we find Note that H depends on (t, y) through the traveling wave coordinate w.
Second, suppose instead that all filaments have the same length ℓ with probability one (see Fig 2e). In this case the only genuine random variables are X 0 and V. Then we have Hence Eq (4) represents the density of fluorescent filaments at time t and location y, or the theoretical profile curves, when the velocity is fixed, and Eq (5) represents the density of the fluorescent filaments (or profile curves) when the length is fixed. The details of these calculations are found in S1 Appendix.

Initial setup
As will be shown in Section 3.3, the mathematical theory indicates that data from FRAP experiments reveals little information about filament length distributions, however some information can be obtained. By knowing how the density of the filaments changes in time and space some limited information about filament length can be deduced. In order to explain this we consider the initial setup for the system and distinguish between the initial distribution of the right endpoints of filaments and the distribution of filament densities. The first is independent of filament length ℓ and the second is not. For our mathematical setup we consider the filaments where the right endpoints are uniformly distributed in the interval [0, F] (see Fig 1a). In addition, if the filaments are long, the density measured from data will be smoother (less variability, with fewer filaments) but if the filaments are short, the measured density will be noisy (more variability which will require more filament measurements to smooth the density signal) as seen in panel (b).
To summarize, filaments with a short average length will give density measurements with greater variability and spatial derivatives. Measurements near the cell membrane may give an idea of average filament length due to the boundary.

Type 1 simulations and theory
In type 1 simulations the velocity for each filament does not change with time (see Fig 2). First we consider a special case of this, namely when the velocity is the same for all the filaments and the filament length, which does not change during the course of the experiment (here we do not model filament assembly or disassembly), is random. Then we consider fixed lengths and random (but constant) velocities. Finally we consider both random velocities and random lengths.
3.3.1 All filaments have the same constant velocity. For fixed velocity, Eq (4) shows that the profile curves are traveling waves since they depend only on the traveling wave coordinate w = y − vt. In other words, the shape of the profile curve does not change but rather is shifted to the right over time. When the bleached region is where the filament density is constant (see Fig 3a and 3b) the profile of the traveling wave will not give much information about the length distribution. Recall that we assume the right endpoints of the filaments are uniformly distributed in the region [0, F]. Because we are placing the right endpoints in [0, F], the density of filaments ramps up from zero to a constant, remains constant from zero to some value below F, and then ramps down to zero at F (see Fig 3a and 3b). Thus Eq (4a) relates to the region left of the right endpoint placement (for instance Fig 3a left of 0) where the density is zero or ramping up. Eq (4b) describes the wave in the rest of the domain up to the bleached region, and Eq (4c) describes the wave to the right of the left edge of the bleached region. We consider only positive velocities and values in the bleached region, the main variations of concern are the transition from the non-bleached region to the bleached region, that is values near y 0 . For values of w such that F − w > u where u is a value such that R u 0 dm L � 1, the first integral in Eq (4b) dominates and is the mean filament length. In this scenario, the wave front is far enough away from F so the profile curve is constant (the boundary effects are negligible-see Figs 3 and 4). Thus the profile curve will be a wave which jumps down from a constant value (determined by the first integral in Eq (4b)) to zero (when Eq (4c) is used) and the only information about the length distribution that can be determined is the integral condition given above which says something about the length of the interval "containing" most of the density of L.
As an example see the first row of Fig 4. In Fig 4a, curves of a scaled version of Eq (4) are shown with data from simulations where the bleached region is in the region where the filament density is uniformly distributed (not just the right endpoints). The filament density is shown for the full domain before bleaching in the inset. The region to be bleached, between y 0 = 50 and y 1 = 60, is shown as grey. The profile curves are constant with a jump at the transition to the bleached region. In Fig 4b the  To summarize, typically, if the filament velocities are constant the wave profile will have a front at y 0 which moves forward into the bleached zone without changing shape and there is almost no information about the length distribution. It may be possible to learn something about the filament lengths if the bleach region is near the cell membrane where the filament density may not be constant due to the boundary imposed by the membrane.
The fluorescence intensity profiles obtained experimentally are not traveling waves (see Fig  1c and 1d). They have an abrupt transition at the time of bleaching from the fluorescent region to the bleached region, but as time evolves the transition becomes smoother and less abrupt. In order to explain the experimental data, we explored the effects of random filament velocity on the curves. From now on, the bleached region will always be in the "plateau" region where both the right endpoints of filaments and the filament density are uniformly distributed to avoid boundary effect in the theoretical and simulated results. Since we only consider filaments moving to the cell periphery, the right side of the profile curves and simulations do not give any additional information.
The rest of section 3.3 will show a primary result of this work-that the fluorescence intensity profile is affected only by filament velocity, and not by filament length.

Filaments with random velocity and fixed length.
The filament velocity distribution affects the fluorescence intensity profiles profoundly. The equation derived for fluorescent filament density with random velocities and fixed length filaments (Eq (5)) is compared to numerical simulations of type 1 in Fig 4d-4l. When the velocity is uniformly distributed the profiles are piecewise linear with the slope changing with time. For Gaussian distributed velocities the profiles are smooth curves with more abrupt transitions than the uniformly distributed velocities. The profile curves are also smooth curves in the case of gamma distributed velocity but the transitions are not as abrupt as when the Gaussian distribution is used.
We then turned our attention to how changes in the fixed length parameter affects the results of type 1 simulations. It is clear from Fig 4d-4l that as the length of the filaments increases the profiles from simulations approach the theoretical curve which is determined by the velocity distribution. If the length is small, there is more variation in the simulated results due to the random nature of the simulations but that variation can be averaged out giving profiles which are similar to those with filaments of longer lengths. Fig 5c-5e, shows how the standard deviation of type 1 simulations vary with the filament lengths. In these simulations (similar to those shown in Fig 4) the averages of 50 simulations are plotted with error bars indicating the standard deviation for simulations with filaments of varying lengths and with uniformly distributed velocities. Clearly the variation increases as the filament length decreases. Depending on the quality of the data it may be able to surmise length information based on the noise in the data.

Filaments with random velocity and random length.
Next, we allowed the length to vary according to different distributions with the same mean while keeping the velocity distribution fixed as a gamma distribution. We do not have a theoretical curve to compare with these type 1 simulations because both velocity and length are random. Fig 5a, curves obtained with the three different distributions (uniform, Gaussian and gamma with the same mean length) for lengths are almost superimposed. In panels (f)-(h) the length distribution was left the same (a uniform distribution) but the mean length was altered from 0.5, 2, 4 microns. Again when the filaments have shorter lengths the data is noisier but still follows the same basic curve (determined by the velocity distribution).

Comparing type 1, 2, and 3 simulations
Here we compare results from type 1, 2, and 3 simulations with fixed filament length and various random velocities. Recall we simulated the velocity in three ways: type 1 each filament gets a different velocity but it remains constant throughout the simulation, type 2 velocities change after T on time units have elapsed but always come from the same distribution, and type 3 filaments have velocities which change after T on time units but then they stop in between velocity changes for T off time units.
Simulations of type 1 match the results from Eq (5) (Fig 6a, 6d and 6h). For type 2 simulations, the initial part of the wave front is slower (the filaments with larger velocities on average do not maintain the large velocity and thus they do not move as far into the bleached region as before) but the back of the wave front is faster (on average the filaments with slow velocities do not remain slow). The overall effect of the change is to make the transition from the bleached region to the unbleached region more abrupt than before (Fig 6b, 6e and 6i). This has the largest effect in the case of the uniform distribution where the transition is linear for type 1 but nonlinear and sharper for type 2 (Fig 6b). For type 3 simulations, filaments paused for a period of time before changing velocity. The pausing time follows a uniform distribution with mean 0.5 minutes. The resulting profile curves are similar to the profiles of simulations of type 2 (without the pausing) except the velocity of the profiles is multiplied by the fraction of filaments that are moving, t on t on þt off , where τ on is the mean time the filaments have a fixed velocity and τ off is the mean time the filaments are paused before changing velocity (Fig 6c, 6f and 6j). The profiles are somewhat advanced at the front end and somewhat delayed at the back end. The pausing seems to slightly ameliorate the sharpening of the wave caused by the velocity changes. Type 2 and type 3 simulations give similar results. When filament velocity is a random variable, the simulated profile curves are not traveling waves and the abrupt change at the time of bleaching is smoothed out as time advances as is seen in the experiments.

Fitting theory and simulations to the data
We have a theoretical formula (Eq 5), and three types of simulations to compare to experimental data. We do not compare type 2 simulations to the data for two reasons. First, the results from the last section show that simulations of type 2 and type 3 are similar. Second, for long τ on (i.e., longer than the duration of the experiment) type 2 and type 3 simulations are the same. We now optimize the parameters in the theory, type 1, and type 3 simulations to fit the experimental data.
We fit thirteen data sets from five different experiments. All the data was from polarized migrating astrocytes showing asymmetric profiles of fluorescence recovery indicative of the polarization of the microtubule driven transport of IFs [12].
In all the simulations, filament length is uniformly distributed on the interval [.05, 10] with mean 5.025 microns. For computational convenience we shifted the distribution away from 0. For the theory and type 1 simulations we optimized over the average velocity, μ, for the uniform distribution, the average velocity, μ, and standard deviation, σ, for the Gaussian and the gamma distributions which is parameterized by the shape parameter k and the scale parameter  ) and the yellow x's show results from type 2 for comparison (with the mean velocity 2 3 of the comparable type 3 simulation). In (a)-(c) the velocity is uniformly distributed, in (d)-(f) it has a Gaussian distribution, and in (h)-(j) it has a gamma distribution. All simulations have filaments with length 10 microns. The other parameters are the same as in Fig 4. The y axis is fluorescence intensity (a.u.) in all panels. Length: All-fixed.
https://doi.org/10.1371/journal.pcbi.1010573.g006 type 3 simulations have three to five free parameters. The type 1 and type 3 simulations are stochastic processes so we optimized 30 realizations for each data set. The objective function minimized was where s(t i , x j ) is the fluorescence intensity of the simulated data, d(t i , x j ) is the experimental data, N is the number of time points, and M is the number of spatial points. Fig 7 shows the average of the minimum of the objective function E (best fit simulation) for each of the data sets. Our model comparisons are only based on how well the model data fits the experimental data. As the models considered are not nested, statistical tests are not applicable for comparison. Although, both type 1 and type 3 simulations fit the data well, the best fits come from type 3 simulations as can be seen, for one data set, in Fig 8. Furthermore, in all cases we successfully simulated most of the fluorescence curves using a truncated gamma distribution for the velocities. This is not surprising since the videos of FRAP experiments clearly show a strong disparity in the filament speeds, with only a few filaments moving very fast and a large majority moving very slowly (see videos [12]). The gamma distribution is the only asymmetric distribution considered thus allowing for a fat tail. Using the scenario (type 3 simulations with gamma distributed velocities) representing the best experimental intensity profile curves, information about the filament dynamics can be extracted from the 13 data sets considered. Fig 9a summarizes the filament velocity distributions, the mean velocities of all filaments and only moving filaments (Fig 9b), the filament mean off and on times, and the percentage of stationary filaments (Fig 9c). The mean velocities of moving filaments are found to range from 0.004 to 0.05 microns per second (with an average over the 13 data sets of 0.0194) and the average velocities of all the filaments (including moving and stopped filaments) and range from 0.003 to 0.025 microns per second (with an average of 0.0108). On average the percentage of stopped fibers at steady state is t off t on þt off . Thus the data shows that in nine of the 13 data sets over half the filaments at any time are stationary (Fig 9c).

Discussion
The work here builds a theoretical framework to quantitatively analyse the directed transport of anisotropic structures and allows the reconstruction of fluorescent profile curves generated by FRAP experiments. Here the focus is on IFs but the work is more generally applicable to all kind of anisotropic structures, such as mitochondria which are also actively transported along microtubules [23]. We show that data from FRAP experiments on IF, namely fluorescence intensity profile curves, reveal important information for determining the velocity of the filaments including mean velocity and shape of the distribution. In fact, the most important distribution governing the shape of the intensity profile curves is the distribution of the filament velocity, μ V . Profile curves which are observed to be traveling waves would suggest that all filaments are moving with the same constant velocity. On the contrary, profile curves which are t on þt off . Panel (c) shows the mean on and off times predicted by type 3 simulations for the 13 data sets considered. The squares to the left are the mean off times, τ off , and the squares to the right are the mean on times, τ on . The black circles are the averages of mean off and on times over the 13 data sets. The grey region shows the percent of filaments which are stopped (the horizontal coordinate is randomly perturbed for viewing purposes). Both on and off times are uniformly distributed. The lines connect values from the same data sets. The colors indicate the same data set. Length: All-uniform with mean 5.025 microns.
https://doi.org/10.1371/journal.pcbi.1010573.g009 piecewise linear with the slope changing in time suggest that the filament velocity is constant in time and uniformly distributed. Profile curves which are sigmoidal indicate velocities are normally or gamma distributed. Finally, we found that the filament length distribution has no impact on the global dynamics of filaments.
Our results show that for polarized migrating astrocytes with reduced retrograde transport, a gamma distribution for the velocity of the filaments best matches the data. Allowing the filaments to pause and restart with new velocities, gives the best experimental fit. In fact, in most of the data sets over half the filaments are stationary and most of the moving filaments have a small velocity. However, some filaments have a large velocity as indicated by the long tails in the velocity probability densities in Fig 9a and the long tails in the profile curves in Fig 8. This is consistent with IF transport by one directional motor molecules with friction and with IF experiments where velocities were described to be non Gaussian with a high propensity of slow filaments in previous work [19,20,24]. Moreover, average velocities extracted from the fit of the profile curves with gamma distributions are consistent with values published in the literature [25], but lower than the range described for the transport of isolated small filaments called squiggles [25,26] and longer filaments measured in [17]. The distributions are consistent with filaments having a large range of velocities which is seen experimentally [25]. All of these velocities are lower than the maximum velocity of an individual kinesin-1 molecule of 0.668 microns per second [27] but which depends on the concentration of ATP and can be as low as 0.045 microns per second for low concentrations. Additionally, increasing the load on the kinesin molecule will decrease the velocity [27]. Long pauses have also been observed to last up to 25-80 percent of the time [25]. This matches our model predictions that at any one time fewer than half of the filaments are moving and the average percentage of time the filaments are paused ranges from 25 to 73 (with an average of 53.5).
The very slow speeds are not surprising because there are a lot of sources of friction for the transport of IF mainly driven by kinesin motors: interaction with other organelles, crosslinking proteins between microtubules and intermediate filaments such as plectin, or dyneins, although dynein activity has been shown to be inhibited [12]. Maintaining a network of filaments that is constantly being restructured requires a delicate balance where a portion of the network is stable, a portion is being dismantled, and a portion is being constructed.
Our results in the case of one directional transport also predict that other types of profile curves can be recovered with velocities which are not gamma distributed. These predictions could be used in other types of experiments investigating different cell conditions not described in [12]. For example, when the friction in the system is reduced, we expect the velocity distributions to be more symmetrically distributed (Gaussian) because only one motor is involved. This could be achieved by reducing the cross-linking of the filaments. Two possible methods would be to inhibit plectin, or completely inhibit dynein. Uniform and Dirac velocity distributions are ideal cases which were examined for the sake of comparison and are less biologically motivated.
There are several possible causes that would inhibit filament motion. Physical obstacles could hinder filament transport including the crowding from other filaments or crosslinking to other filaments. Stalled velocity due to the tug-of-war caused by motor molecules of opposing types could be another possible reason. Mathematical modeling shows that there are several scenarios where the majority of filaments remain in a state where the tug-of-war is unresolved resulting in stalled filaments [15]. Filaments detaching and staying detached from the microtubule would also be stalled. Finally there are direct and indirect interactions with actin which could cause anchoring of the filaments [17]. Overall, the data suggests that many of the fluorescent filaments are in the static portion of the network or in the process of being disassociated or associated with it. Thus the majority of filaments are stationary or moving with very slow velocities. This does not preclude the possibility that many filaments which are not associated with the stationary portion of the network are also pausing for long periods of time.
Assembly and disassembly of IF occurs on a time scale of hours in neurons [28]; whereas, the FRAP experiments take place on the time scale of minutes. Yet in epithelial cells, keratin assembly/disassembly occurs at the time scale of minutes [29]. But in our experiments the profile curves remain low in the middle of the bleached region indicating that assembly is not playing an important role. Preliminary results of simulations with length changes suggest that the results presented here are robust. This is not surprising since length does not substantially affect the profile curves. Of course, in systems where filament length can change at a time scale comparable to that of active transport, our analysis is not appropriate and the profile curves will have different characteristics.
Let's consider three possible refinements to model. First, we could allow filament transport in two directions. This would be important when considering symmetric profile curves, for example those observed in astrocytes 8 hours after wounding, when cell polarity is well established [12]. Two directional analysis would indicate whether the filament velocity has the same characteristics in both directions. When there are no non-fluorescent moving filaments, the bleached region is symmetric, and the velocities are equal in magnitude but opposite in direction no new information is gained. If there are non-fluorescent, moving filaments the profile curves will be distorted but in a symmetric manner. If the velocity distributions are the same but with different mean velocities (still in opposite directions) the moving non-fluorescent filaments could break the symmetry. Second, we could allow the filament velocity to be correlated with filament length [14]. How this correlation would alter the profile curves is hard to predict without knowing how the length and speed are related. Third, we could consider the elastic nature of IFs. The elasticity of the filaments has two possible relevant effects: length change and speed change [14]. The first should not affect the profile curves but the second could. Finally, by combining the information about filament velocities found here with models of filament transport [14,15], it should be possible to elucidate properties of motors involved in the transport. Additionally, when the cell is stationary, while the filament network is very dynamic, there is no net change in filament density.
In future mathematical models, we plan to investigate conditions which would allow a steady state density and give insight into how the cell maintains a dynamic network with no net filament transport. Further mathematical analysis could be done using other mathematical formalization. A standard method would be to use the Chapman-Kolmogorov equation, but for our problem that requires simplification. All three types of simulations can naturally be thought of as realizations of stochastic processes. With some work and simplifying assumptions, they can be framed as Markov processes. Type 1 simulations, while certainly stochastic, have all of the randomness front-loaded at time 0. Consequently, they are trivially Markovian because the state of the system at each strictly positive time depends deterministically on the state at any previous time. The easiest way to frame type 2 and type 3 simulations as Markovian is to consider only the velocity process and to require the elapsed time between velocity jumps to be exponentially-distributed (rather than uniform), leaving us with classical examples of continuous-time, time-homogeneous Markov processes. For interested readers the resulting Chapman-Kolmogorov equation and other related formulas are given in S2 Appendix.
In summary, the modeling framework proposed in this work provides an in-silico platform to study the impact of IF protein post-translational modifications [30] or mutations [31], depleting or silencing one motor type, or altering the IF network composition on the IF transport and organization in cells.