A DRIFT-DIFFUSION MODEL FOR MOLECULAR MOTOR TRANSPORT IN ANISOTROPIC FILAMENT BUNDLES

In this study we consider the density of motor proteins in filament bundles with polarity graded in space. We start with a microscopic model that includes information on motor binding site positions along specific filaments and on their polarities. We assume that filament length is small compared to the characteristic length scale of the bundle polarity pattern. This leads to a separation of scales between molecular motor movement within the bundle and along single fibers which we exploit to derive a drift-diffusion equation as a first order perturbation equation. The resulting drift-diffusion model reveals that drift dominates in unidirectional bundles while diffusion dominates in isotropic bundles. In general, however, those two modes of transport are balanced according to the polarity and thickness of the filament bundle. The model makes testable predictions on the dependence of the molecular motor density on filament density and polarity.

1. Introduction.Cells depend on intracellular transport of cargo like vesicles, organelles and protein complexes.Simple diffusion is rarely effective; usually, this transport is driven by the movement of processive, cargo binding molecular motor proteins along polar filaments, such as microtubules and actin [5].Examples are motor proteins of the dynein and kinesin families which move along polar microtubule fibers to minus and plus ends, respectively [15] (figure 1), and proteins of the myosin family transporting cargo along actin filaments [8].
The polarity, i.e. orientation of plus and minus ends, of filaments within the bundle in general changes along the bundle ( [14]).Our focus is the effect of variations of polarity and bundle width in space on the distribution and dynamics of molecular motors.It is known that molecular motors frequently detach from one filament and attach to another, and as a result frequently change their direction of movement [9].As a result, their density therefore would evolve according to a diffusion-drift process rather than through mere drift.
While this case has not been treated in the modeling literature, a drift-diffusion model was formulated in [13] that relied on the assumption the filaments are long, spanning the whole filament bundle, while diffusion is caused by random motion within the cytoplasm.A variant of this model has been analyzed using variational methods [3] and various other variants of it have been used to describe unidirectional transport in axons coupled to backward diffusive flux within the cytoplasm ( [11,10]).Recently a comprehensive review on modes of intracellular transport ( [2]) has been published which thoroughly elaborates on the interplay of directed transport along filaments and diffusive modes of transport in the cytoplasm.We also note that the deep connection between transport models and drift-diffusion equation was elucidated by Peter Lax and colleagues [4], and the general idea of deriving a drift-diffusion model from more complex transport models has a long and venerable tradition in mathematical biology [7,6,12].
Our aim here is to derive the drift-diffusion equation for molecular motors moving along anisotropic filament bundles.We are considering generic molecular motor proteins which move along protein fibers at a given velocity.The fibers themselves are assumed to have the same characteristic length and to be aligned in parallel forming a bundle that corresponds to a microtubule-filled axon [1] or to a bundle of actin filaments coated with myosin motors engaged in cargo transport [17].The filaments are assumed to be polar and might therefore, within the filament bundle, point in one of two possible directions (figure 1).We call those two fiber ends minus end and plus end and assume that our generic molecular motors move towards the plus end at a given fixed speed.
Figure 1.A: Cargo bound molecular motor proteins move along filaments.Arrowheads mark the polarity of the filaments which determines the direction of movement.B: Once molecular motors reach the tip of the filament they, a), fall off and, b), immediately reattach to a another filament that is randomly picked out of those, which overlap with the molecular motor protein.
We further assume that the motors are processive and move to the filament ends; once a molecular motor reaches the plus end, it falls off that specific filament and immediately reattaches to another filament that is picked randomly out of those which overlap with the instantaneous molecular motor position (see figure 1).As that randomly picked filament might have any of the two polarities, the molecular motor will either keep its direction or turn and move in the opposite direction.
It is clear that molecular motors within a bundle will undergo frequent pauses and/or changes of direction, and their distribution will evolve according to a general drift-diffusive pattern.If one filament polarity dominates, we should expect more drift; if average polarity is zero, we expect more diffusion.One could therefore come up with an ad hoc model that combines both effects, weighting them with coefficients that depend on the anisotropy of the filament bundle.
In this study it is our aim to derive such macroscopic drift-diffusion model as perturbation approximation of a detailed transport model which includes information on the position of molecular motor binding sites on specific filaments and their polarities.To this end we introduce a scaling of respective parameter with regard to the filament length and we derive the first order perturbation model as an approximation of the original microscopic model in the case of slowly varying polarity and width of the filament bundle.
The result is a drift-diffusion model equation which includes expressions for drift and diffusion weighted by local bundle polarity, but it also includes additional velocity fields which are caused by changes in bundle width and polarity, which would otherwise be overlooked in an ad hoc model.The resulting model for the density of motor proteins ρ = ρ(t, x) where x ∈ R corresponds to space and t ≥ 0 to time reads where V m > 0 is the constant speed at which molecular motors move along filaments and l > 0 is the length of filaments.Here we omit initial and boundary data since those typically need to be adapted to the specific modeling application.The distribution of plus-end-forward and minus-end-forward directed fibers is given by r + (x) and r − (x) respectively.Their sum r(x) = r + (x) + r − (x) represents the width of the filament bundle, i.e. the number of filaments per cross-section, while their difference r(x) = r + (x) − r − (x) represents the polarity of the bundle.The relative polarity ranges between −1 and 1 and the order parameter Q 2 represents the anisotropy of the bundle and ranges between 0 (isotropic) and 1 (fully anisotropic).
In addition, we also take into account the case where the two types of oppositely oriented filaments with densities r + = r + (t, x) and r − = r − (t, x) are moved by given velocity fields v + = v + (t, x) and v − = v − (t, x).Notation-wise we write their average velocity as v = (v + + v − )/2 and the relative velocity as v = (v + − v − )/2.The modified evolution equation for the density of molecular motors in this case reads The average velocity affects the dynamics through an additional drift term while the relative velocity v acts as a correction to the molecular motor speed V m .
The structure of this papers is as follows.In section 2 we formulate the microscopic structured model which describes the population of molecular motors retaining the information about the polarity of the filament to which they are attached and the relative positions of the binding sites.
In section 3 we derive a first order perturbation equation which represents an approximation to the original model under the scaling assumption that the filament length is small compared to the bundle length scale.
In section 4 we analyze the relative magnitude of the diffusive and the drift terms in the resulting drift-diffusion model and demonstrate the effect of one of the drift terms using numerical steady state solutions.Note that as the external velocity fields only add correction terms we restrict the final discussion to the equation for a static bundle (1).Finally section 5 will be devoted to concluding remarks.Figure 2. A: A molecular motor at position x moves at speed V m in the direction of the filament plus end marked with an arrowhead.Its relative position η ∈ [−l/2, l/2] represents the distance towards the filament center point c, hence η = x− c. η = l/2 corresponds to the molecular motor being positioned at the plus end and η = −l/2 refers to the minus end.B: The same situation for a molecular motor attached to a minus-end-forward direction filament.The relative position is now computed as c−x and therefore the filament center is at x + η while η = l/2 and η = −l/2 still refer to the plus end and minus end, respectively.
2. Microscopic model formulation.First, we formulate a structured reactiontransport model which describes two types of molecular motors, bound to plusend-forward and minus-end-forward filaments, respectively.The motor densities χ ± = χ ± (t, x, η) are interpreted as densities with respect to the position within the filament bundle x and with respect to the microscopic dimension η ∈ [−l/2, l/2] which corresponds to the distance of the molecular motors to the center point of the filament they are attached to.The interpretation of η varies between the two groups of molecular motors.However, η = l/2 always refers to the plus end of the filament and η = −l/2 refers to the minus end in both cases (see figure 2).
The equation for χ The transport term in the direction of η models the fact that molecular motors move towards the plus end with speed V m .Transport in physical space is additionally affected by transport of the fibers themselves and therefore the speed with respect to physical space is v + (t, x − η) + V m where x − η is the position of the center point of the filament while the molecular motor itself is at position x (see figure 2).Molecular motors detach from the fiber at η = l/2, i.e. once they reach the plus end.The total number of molecular motors detaching from their filaments per unit time at x is given by V m (χ + (t, x, l/2) + χ − (t, x, l/2)).The molecular motors immediately leave the minus end of filaments, therefore the boundary condition at η = −l/2 is zero density.Detached molecular motors reattach with equal probability to any filament within the bundle cross-section at x.This includes actin filaments with center points between x−l/2 and x+l/2 and therefore molecular motors attach to filaments at any relative position between η = −l/2 and η = l/2.The polarities of filaments, however, vary according to the concentrations r + and r − .Therefore reattachment to filaments of positive polarity at x − η, or of negative polarity at position x + η is weighted the coefficients R + and R − respectively, which have the form Note that R + and R − are probability densities in the sense that l/2 −l/2 (R + +R − ) dη = 1 and therefore reattachment conserves the total number of molecular motors.
The analogous reaction-transport model for the molecular motors moving in negative direction is given by where the effect of the molecular motor speed V m has opposite sign in the transport term with respect to physical space but retains the same sign in the transport term with respect to the position relative to the filament center.
With a view to derive one single equation for the combined density of rightmoving and left-moving molecular motors we recast this system of equations for the new quantities χ = χ + + χ − and χ = χ + − χ − , which represent the total density of molecular motors and their polarity.We also recast the velocity fields according to and we obtain the rewritten system of equations taking the sum and difference, respectively, of the system (4), (6), 3. Derivation of the first order perturbation solution.Based on the system of equations ( 7), ( 8) we consider the asymptotic regime where the filament length l is small compared to the length scale on which the bundle polarity and width along the bundle vary, and we derive the first order perturbation equation as an approximating model.
3.1.Scaling.We scale the distance between the molecular motors and the center of filaments by the length of the filaments, η = η/l and we drop the tilde to keep the notation simple.Note that now the boundary value of χ is given at η = −1/2 and the domain of the microscopic dimension is the interval [−1/2, 1/2].The characteristic scales for x and t are ∆x (on which the polarity and width of the bundle change) and t 0 ∼ ∆x/V m .What we have in mind is the asymptotic regime where l = l/∆x → 0. Note that implicitly we also scale the densities χ = χ 1 l∆x .First, we scale the probability densities for reattachment (5), replace r ± (x ∓ lη) by Taylor expansions at x and thus find asymptotic expressions for the sum and difference of R + and R − which appear in the system ( 7), (8), Next, we scale the system (7), (8), plug in the asymptotic expansions (9) and also replace v(x ± lη) and v(x ± lη) by Taylor expansions, which leads to Finally, we introduce the asymptotic expansions of χ and χ and for the macroscopic densities with respect to x we write Plugging in ( 12) into ( 10) and (11), we obtain equations of various orders with respect to l, which we use in the sequel to assemble an O(l) perturbation equation.

3.2.
System of order l −1 .The equation of order l −1 contains the information on how χ 0 depends on η, It implies that and therefore we obtain, by integration on [−1/2, 1/2], the expression Likewise the equation of order l −1 for χ0 reads based on which we conclude by integration that μ0 = 1 2 r r χ 0 (t, x, where we used (14).

3.3.
Equations of order l 0 .The order l 0 equation of (10) reads It contains information on the macroscopic dynamics, namely the 0-order perturbation problem, which we obtain directly by integration and using ( 16), To isolate the microscopic information contained in (17) note that due to ( 14) and ( 16) it holds that and therefore (17) implies that which allows to infer by integration that The order l 0 equation of (11) on the other hand reads The macroscopic part of this equation, due to ( 13) and ( 15), can be written as for a function C(t, x) that does not depend on η and for which we find by integration that For the profile of χ1 with respect to η this means according to (20) that and therefore we conclude Finally integration on [−1/2, 1/2] yields an expression for μ1 which corresponds to the order-1 version of ( 16) where we used previously gathered information, namely (19), (21) as well as ( 14) and ( 16).
3.4.Equation of order l.Summarizing the coefficients of order l in (10) we obtain (23) After integration with respect to η the right hand side cancels and we obtain the evolution equation for µ 1 , where we might substitute (24) to eliminate any dependence on the "bar" quantities μ0 and μ1 .
3.5.First order approximating problem.We multiply (24) by l, add the 0order limit equation (18) and substitute (22) obtaining We aim to convert (25) into a first order perturbation equation for ( 10), ( 11) by adding appropriate expressions of order l 2 .First, we replace any occurrence of µ 0 + lµ 1 in (25) by ρ.We also replace any occurrence of µ 0 alone by ρ which implicitly means that we added the missing multiples µ 1 which are of order l 2 .Thus we obtain Note that if we substituted the expansion ρ = ρ 0 + lρ 1 + ... into (26), the equations for ρ 0 and ρ 1 would coincide with (18) and ( 24) respectively and therefore the solution to (26) is a first order approximation to µ = 1/2 −1/2 χ dη.The mixed derivative ∂ xt ρ in (26), however, renders any straightforward interpretation of that equation difficult.We therefore undertake a second step in which we again add terms of order l 2 which will allow to rewrite (26) as a classical driftdiffusion equation.Note that the structure of (26) corresponds to where U 0 , U 1 and G do not depend on ρ.The approach which we take is to add the second order expression −l 2 ∂ x (G∂ x (G∂ x (U 0 ρ))) which allows to write (27) as The rewritten equation ( 28), however, is the sum of one equation and its derivative with respect to x multiplied by l.Hence the solution of the following simpler equation 0 ) also solves (28) and coincides with the original equation ( 27) up to the first order terms.
This procedure applied to (26) yields where we used the notation ( 2) and where we already distributed expressions to the right and to the left of the equality sign according to their diffusive, respectively drift-type nature.Note that setting the fiber velocities v and v to zero we already recover the first order perturbation drift-diffusion equation for molecular motors in a static filament bundle (1).Lastly, we will derive the evolution equation for Q and substitute for ∂ t Q in (30) which will eliminate many of the expressions in (30) which depend on v and v.We use the fact that the velocity fields v + and v − transport the filament fibers, hence it holds that which in v and v notation and using r = r + and r = r + − r − is equivalent to From this we derive which we plug into (30) to obtain the first order perturbation equation 4. Discussion.The first order approximation model (31) has the classical structure of a drift-diffusion equation.For the further discussion we set the velocity fields of the filaments to zero as they only introduce correction terms.For a static bundle the model reads (32) where we introduced abbreviating notation for the velocity fields and the diffusion coefficient.This equation contains the dominant drift term with corresponding velocity V 1 and the diffusive term with corresponding diffusion coefficient D. There are also small drift terms with corresponding velocities V 2 and V 3 .All three drift terms and the diffusion term depend on the local polarity Q(x) of the bundle.
First, let us analyze the relative magnitudes of the three drift terms and the diffusion term.Recall that ∆x is the characteristic length on which polarity Q(x) and bundle width r(x) change, and that the scaling assumption under which we derived (32) was l ≪ ∆x.In general it holds that Q ∼ 1, so the magnitudes of the diffusion coefficient and the velocity fields, under conditions which are compatible with the scaling assumption, are .
In an almost unidirectional bundle, where Q 2 ≈ 1, V 2 and D that are proportional to (1 − Q 2 ) are negligible, while V 3 ≪ V 1 for any value of Q, and the model reduces to purely drift-like transport with velocity V m .The other terms start to contribute significantly in an almost isotropic (or, as it is called in the biological literature, bipolar) bundle, when |Q| ≪ 1.Let us first consider the case when r(x) = const.Then, in the zero flux stationary case, the motor density is governed by the equation the solution of which is This solution predicts that the motors will be distributed smoothly in the almost isotropic bundles with maxima and minima of their densities corresponding to zeros of polarity Q(x).If we neglected the diffusion term, the motors would concentrate into delta-function-like aggregates at zeros of polarity Q(x), so the diffusion term has non-trivial smoothing effect.Note that even in a highly anisotropic bundle this term is not negligible in the vicinity of zeros of polarity Q(x).
In situations where the width of the filament bundle r varies, the velocity field V 2 states that mere thickening of the filament bundle triggers a drift in the direction of the thickening.This effect is small in the highly anisotropic bundle when Q 2 ≈ 1, but is not negligible in an almost isotropic bundle when |Q| ≪ 1.The fact itself that thickening of the bundle should cause directed transport of molecular motors is somewhat counter-intuitive, but it can be explained by the following argument.In the example of figure 3 sudden thickening of an isotropic bundle happens within Figure 3. Schematic illustration of bundle thickening causing a net drift of molecular motors: In this example sudden thickening of an isotropic bundle happens within a narrow transition zone.Once molecular motors attach to a filament near to its minus end they leave the transition zone along this filament.The probability to leave the transition zone in the direction of the thicker part of the bundle, however, is higher than in the other direction which translates into a net drift in that direction.a narrow transition zone.Molecular motors might circle back and forth within the transition zone as long as they hop onto filaments near their plus ends.Once they attach to a filament near to its minus end they leave the transition zone along this filament.The probability to leave the transition zone in the direction of the thicker part of the bundle, however, is higher than in the other direction for the simple reason that more filaments enter the transition zone from the side where the bundle is thicker.In the example of figure 3, thickness increases by the factor 2 and therefore also the number of molecular motors leaving the transition zone in the direction of the thicker bundle is larger by that same factor.It is clear that this effect must go away as the bundle becomes more anisotropic and also as the filaments become shorter which is reflected by corresponding factors in V 2 .
In the special case of perfectly isotropic bundle where Q ≡ 0, a stationary solution corresponding to zero motor flux can be found from the equation the solution of which is: ρ(x) = C √ r where constant C > 0 is given by a normalization condition.Hence we predict that in isotropic filament bundles of varying filament density per cross-section, the concentration of molecular motors per unit length is proportional to the square root of the density of filaments per cross-section segment.
The velocity field V 3 represents a net drift in the direction of increasing anisotropy Q 2 and it can therefore accelerate the zero order drift V 1 or slow it down.It should be noted that this effect is restricted to anisotropic bundles due to the factor Q in V 3 .In terms of magnitude, this velocity field is always dominated by the O(1) ≪ 1 and we conclude that V 3 can be dropped from (32).
Although its contribution is minor we show numerical steady state solutions below to illustrate the effect of V 3 in scenarios which, rigorously speaking, violate our scaling assumption l ≪ ∆x.
We consider a filament bundle of constant width and given a specific polarity pattern and use the parameters V m = 0.5 µm/s and l = 30 µm. Figure 4 shows the stationary solution for a given density at the left boundary given by ρ 0 = 0.1 µm −1 and with the polarity Q as shown in the lower picture.Here the molecular motors, since polarity is zero in the left part of the bundle, enter the bundle by diffusion and since the right part of the bundle is unidirectional they are transported out of the bundle at the maximal speed V m .In the center there is a transition zone where the polarity grows linearly.
The first graph in Figure 4 shows the stationary solution of (32) while the second graph shows the stationary solution of the equation after we have dropped V 3 .In that second scenario the stationary density of molecular motors is continuous and monotonic while the V 3 drift term in the first graph accelerates transport within the transition zone.However the stationary density rebounds when the polarity reaches its maximum as the effect of the additional drift isn't effective any more.Note that this rebound effect increases even slightly if the transition zone gets more narrow.However, the effect goes away as the transition zone widens which flattens the polarity gradient and corresponds to the parameter regime which is actually covered by the scaling assumption.
In a second scenario we look at the concentration of molecular motors at the interface where the polarity of a filament bundle suddenly switches sign from positive to negative.We assume that the filament bundle is unidirectional at both ends and that the change of polarity is linear within a narrow transition zone where molecular motors concentrate.In figure 5, the width of the transition zone is 30 µm which corresponds to the filament length.The major difference between the two stationary density profiles with and without the polarity gradient-generated drift V 3 is that the solution of the full equation (32) seems to enforce a flat plateau in the motor density along the transition zone while without that drift term there is a more pronounced peak in the center of the transition zone.Of course, if the transition zone is even more narrow, the stationary concentration profile of the microscopic model would be wider, while the stationary concentration profile predicted by (32) would be supported within the transition zone and therefore correspond to a bad approximation.Thus, when polarity pattern changes on a very short length scale, the drift-diffusion approximation to the original microscopic model is not valid.
Certainly there is a plethora of random effects that causes diffusivity in cargo transport.One is the randomness of observed molecular motor velocities which only very roughly can be assumed to be constant.It was the aim of this study to focus on the diffusive effect of bundle polarity and therefore we made this simplifying assumption.As a matter of fact, the diffusive effect due to filaments switching directions, which is included in our model, is large: the characteristic diffusivity is l/3 V m , where l is the filament length.Comparing this diffusivity to the effect of variations in motor velocity, one finds that, for example, the standard deviation of kinesin velocities is in the range of 10% [16] of the mean motor velocity.Hence, even if motors maintain their individual velocity along their entire runlength, the diffusivity due to variations in velocity between motors will be roughly 10% of the diffusivity effect caused by the filaments switching directions.Similar considerations apply with respect to the diffusive effect of randomness of fiber trajectories.If fibers are transported by molecular motors in a coordinated manner, captured in our model by given velocity fields v + and v − , only deviations from the mean rate of transport would enter the diffusivity.If the single fibers are uncoupled, that effect could be much larger, since motors could transport single fibers in random directions.This situation, however, is definitely not covered by our model.Implicitly, our model assumes significant cross-linking between fibers, which diminishes shear movements of the neighboring parallel fibers, and so the dominant effect is relative sliding of the anti-parallel fibers.

Conclusion.
In this study we derived a drift-diffusion model for the dynamics of molecular motors moving in anisotropic filament bundles.The resulting equation contains three drift terms and a diffusion term that are functions of polarity and thickness of the bundle.In highly anisotropic bundle, the first drift term with the velocity proportional to the bundle polarity dominates.In almost isotropic bundles, diffusion term and drift term depending on the bundle thickness gradient are not negligible and predict smoothing of the motor distribution and motor accumulation in thicker parts of the bundle, respectively.The model makes testable predictions on how motor density per bundle segment scales with the number of filaments per cross-section in isotropic bundles of varying thickness.Finally, the third drift term gives a higher order drift caused by variations of the fiber polarity which is always small in magnitude.

1 QFigure 4 .
Figure 4. Stationary solution (steady transport from left to right) with (top) and without (middle) the term corresponding to polarity gradient-generated drift.Lower graph: given polarity profile.

1 QFigure 5 .
Figure 5. Stationary solution (concentration at the interface of unipolar bundles of opposite polarity) with (top) and without (middle) polarity gradient-generated drift V 3 .The width of the zone featuring a linear transition between the two unipolar bundles corresponds to the filament length.Lower graph: given polarity profile.