A general model for the motion of multivalent cargo interacting with substrates

2.1. The cargo binding model

We posit that the motion of a substrate-bound multivalent cargo originates from the ability of its interaction sites (henceforth referred to as legs) to repeatedly bind/unbind at different positions on the substrate (figure 1). This implies that this process should be able to explain both the diffusive motion of multivalent ligand molecules interacting with receptor-functionalized surfaces [1] and the processive motion of multivalent cellular components interacting with cytoskeletal polymers [5,7,27,28,31]. Differences in cargo transport would then reflect differences in the spatial distribution of leg–substrate interaction rates.

• Download figure
• Open in new tab
• Download PowerPoint

In between binding/unbinding events, the instantaneous position of the cargo will be determined by a balance of forces. For a cargo of length L with N legs, of which n are bound at positions {x_l(t)} at time t, we will define the centre position x_a(t) as the average of the bound leg positions, such that $\sum _{l=1}^n(x_l(t)-x_a(t))=0$. Although this assumption can be generalized, it enforces a simple force equilibrium criterion on the cargo. Each cargo initially binds in the n = 1 state and then unbinds once n = 0, such that its average dwell time while bound to the substrate, t_dwell(x), can be calculated using previously published results [29,30] (see electronic supplementary material, figure S1). Individual legs can bind to (m = 1) or unbind from (m = 2) the substrate with position-dependent rates κ_m(x). However, these rates need to be modified to take into account the size and shape of the cargo. We do this by including two probability distributions S(x|x_a(t), n) and P_b(x|x_a(t), n). The ‘shape factor’ S(x|x_a(t), n) is the probability distribution function that describes the probability of one of the cargo’s N − n unbound legs binding at the position x when the cargo is at x_a(t). This modifies the binding rate to k₁(x|x_a(t), n) = κ₁(x)S(x|x_a(t), n), thereby extending the concept of local concentration of substrate binding sites previously introduced in [33]. For example, a shape factor S_H(x|x_a(t)) = (θ(x − (x_a(t) − L)) − θ(x − (x_a(t) + L))/2L (where θ(x) is the Heaviside step function) states that unbound cargo legs can bind to the substrate uniformly within the range x_a(t) − L ≤ x ≤ x_a(t) + L, defining a putative cargo width of 2L. This approach assumes that legs cannot interact with each other, and that once unbound they can access and try to bind to any position on the substrate allowed by the shape factor, regardless of the time elapsed since their last unbinding. The latter condition requires the typical binding timescale ($t_b\sim 1/max({\kappa }_1(x))$) to be much larger than the timescale associated with the diffusion of unbound legs on the cargo surface (t_d ∼ (2L)²/D_l; D_l is the effective diffusivity of the cargo leg). If needed, these simplifying assumptions can be relaxed and a more complex, possibly time-dependent, shape factor can be introduced to modulate the legs’ distribution around the centre of the cargo x_a(t).

Since cargo legs are stationary once bound, unbinding events can only occur from positions within the set {x_l(t)} (henceforth lengths will be understood to be in units of L). This requirement is simple to implement computationally, but complex to include analytically. Therefore, to make analytical progress, we tackle this problem at the mean-field level by introducing the coarse-grained bound leg distribution P_b(x|x_a(t), n). This probability distribution function describes the probability of one of the cargo’s n bound legs being at position x when the cargo is at x_a(t), and it modifies the unbinding rate distribution such that on average k₂(x|x_a(t), n) = κ₂(x)P_b(x|x_a(t), n). The bound leg distribution is clearly linked to the cargo’s shape factor. For example, if the shape factor only permits binding of cargo legs within a specified range of positions, then the bound leg distribution will decay quickly at the extremities of this range (as will be shown in figure 2a).

• Download figure
• Open in new tab
• Download PowerPoint

In the presence of position-dependent leg binding and unbinding rates κ_1,2(x), a cargo will exhibit position-dependent average displacements and event rates. The ith moment of the cargo displacement distribution, ${\lambda }_m^{(i)}(x_a(t),n)$, and the average rate, ${\overline{k}}_m(x_a(t),n)$, corresponding to each type of event (m = 1, 2) are given by (see electronic supplementary material, §1)

andwhere Δ_m = δ_m,1 − δ_m,2 ≡ 3 − 2m and δ_i,j is the Kronecker delta. The second term in equation (2.1) calculates the mean ith power of the difference between the positions of the binding or unbinding event and the cargo centre, with the factor ∝(1/(n + Δ_m))ⁱ accounting for the smaller displacements of the cargo centre position when more legs are bound. Instead, equation (2.2) is the simple multiplication of the average rate of binding/unbinding over the extent of the cargo by the number of legs available for each type of transition.

We then introduce a third type of transition (m = 3) to take into account any deterministic motion of the cargo independent of position and n. This can be used, for example, to model cargo that bind to substrates that themselves can move [7,28,31], or cargo motion in the co-moving frame of a growing substrate such as the polymerizing tip of a microtubule. This motion can be introduced analytically by assuming that all cargo legs are displaced by Δx (resulting in ${\lambda }_3^{(i)}={(\mathrm{\Delta }x)}^i$) at a rate ${\overline{k}}_3$, and can be generalized to accept different displacement and wait-time distributions.

Armed with these quantities, it is possible to derive the continuity equation satisfied at the mean-field level by P(x, t), the probability distribution function of finding cargo at the position x at time t [34]. A truncated Kramers–Moyal expansion leads to the Fokker–Planck equation (see electronic supplementary material, §2):

where k_on(x) and k_off(x) are the position-dependent effective binding and unbinding rates (respectively) of cargo modelled as single, composite bodies. This effective unbinding rate is equivalent to the reciprocal of the average dwell time of the cargo (see electronic supplementary material, §3). The position-dependent effective velocity and diffusivity, v_eff(x) and D_eff(x), are defined asandwhere P_n(x) is the probability of cargo having n legs bound at position x (see electronic supplementary material, §3; not to be confused with the bound leg distribution), and M = 3 for a system with binding/unbinding events and a net velocity. Together, equations (2.3)–(2.5) predict that cargo can exhibit both diffusive and processive motion depending on the forms of the cargo–substrate interaction rates κ_m(x).

The terms $S_{\lambda }^{(1,2)}(x)$ defined in equations (2.4) and (2.5) both consist of sums over independent contributions from M different types of events, averaged over the expected number of bound legs at each position. They can be understood as the product of a position-dependent total event rate and a suitably weighted average of the displacement moments ${\lambda }_m^{(1,2)}$ (see electronic supplementary material, §2). Substituting equations (2.1) and (2.2) into equation (2.4), we see that the direction of $S_{\lambda }^{(1)}(x)$ depends on $\pm [{\int }_{x_a}|x-x_a|k_m(x|x_a,n)\hspace{0.17em}\mathrm{d}x-{\int }^{x_a}|x-x_a|k_m(x|x_a,n)\hspace{0.17em}\mathrm{d}x]$, which describes the difference in the binding/unbinding rate distributions (± respectively) averaged over the two halves of the cargo width before and after the midpoint x_a. For slowly varying binding/unbinding rates, this has the appealing intuitive interpretation that $S_{\lambda }^{(1)}(x)$ acts in the direction of increasing binding rate (∝k₁′(x|x_a(t), n), where ′ indicates a derivative with respect to x) and decreasing unbinding rate (∝−k₂′(x|x_a(t), n); see §2.3).

This coarse-grained description of cargo motion at the population level can be accompanied by a Langevin equation for the stochastic trajectories x_a(t) of individual cargo (see electronic supplementary material, §5):

Here, dW(t) is the standard Wiener process [35] and the multiplicative noise is intended in the Itô sense. The first term in equation (2.6) defines the effective deterministic motion of a bound cargo due to gradients in the binding and unbinding rates of its legs, while the second term describes the noise affecting its motion. This means that the fixed points of the dynamics will occur at positions x* where $S_{\lambda }^{(1)}(x^{\ast })=0$, meaning thatFor example, for EB-mediated cargo transport in the rest frame of a microtubule growing with velocity v_MT, $-{\overline{k}}_3\mathrm{\Delta }x=v_{\mathrm{MT}}$, and a fixed point will then arise when the effective velocity generated by the binding/unbinding dynamics of cargo legs—the left-hand side of equation (2.7)—is equal to v_MT. Fixed points will be stable when ${(\mathrm{\partial }{S}_{\lambda }^{(1)}(x)/\mathrm{\partial }x)|}_{x=x^{\ast }}<0$, and will then correspond to the positions where cargo accumulate. In the previous example, this would correspond to the cargo tracking the polymerizing microtubule tip.

2.2. Probing cargo dynamics

Having defined the model, we can start probing cargo dynamics by comparing the more intuitive analytical solution with a numerical one. For simplicity, in the following sections, we will assume the uniform shape factor S_H(x|x_a(t)). The bound leg distributions P_b(x|x_a(t), n) cannot be derived analytically in general, but they are straightforward to obtain from stochastic cargo binding simulations. Example bound leg distributions for use in equations (2.1) and (2.2) are shown in figure 2a and electronic supplementary material, figure S4a, together with a heuristic fit which provides a parametric family of analytical curves that can be used in the integrals (see electronic supplementary material, §6). For uniform binding/unbinding rates (κ₁ = κ₂) and ${\overline{k}}_3=0$, figure 2a shows that bound leg distributions are symmetric, and equations (2.1), (2.4) and (2.5) immediately imply that ${\lambda }_m^{(1)}(x_a(t),n)=v_{\mathrm{eff}}(x)=0$. The cargo should then exhibit diffusive motion with no net drift, in agreement with stochastic simulations (see electronic supplementary material, figure S2a).

Interestingly, for N > 2, simulated cargo exhibit distinct short- and long-time diffusivities as a result of their binding dynamics (see figure 2b; electronic supplementary material, figure S2; for 3 ≤ N ≤ 10 each pair of the short- and long-time diffusivities are separated by at least $10max({\zeta }_{\mathrm{short}},{\zeta }_{\mathrm{long}})$, where ζ_short,long are the errors associated with the short- and long-time diffusivities in figure 2b and electronic supplementary material, figure S2, respectively). Cargo always initially bind in the n = 1 state and, for those that survive long enough, the probability P_n(x) of having n legs bound will take some time to relax to the steady-state form (see electronic supplementary material, figure S3). Given that the MSD due to binding/unbinding events depends on the average number of bound legs (for example equation (2.1) shows that ${\lambda }_1^{(2)}(x_a(t),n)\propto 1/{(n+1)}^2$), the diffusivity measured a short time after initial attachment is expected to be higher than that measured later for longer living cargo. Intuitively, the displacement caused by the binding or unbinding of individual legs will be larger when the number of bound legs is small. Figure 2b and electronic supplementary material, figure S2, show that this is indeed observed in simulations, and we propose that it could also explain previously published experimental MSD data showing multivalent cargo diffusivities that decrease as a function of time spent bound [5]. The definition of short and long time is relative to a ‘separation timescale’ $\lesssim t_c=1/{\kappa }_2$ for N > 2 (see the crossover point in electronic supplementary material, figure S2b, and the plateau in electronic supplementary material, figure S3). This separation timescale diverges at N = 2, such that the short- and long-time diffusivities coincide in figure 2b. Figure 2b also shows that the analytical model defined in equations (2.1), (2.2) and (2.5) generates effective diffusivities in agreement with the long-time values obtained from simulations, but only after introducing the steady-state bound leg distributions shown in figure 2a. Excluding the bound leg distributions instead leads to a higher effective diffusivity due to the overestimation of the probability of finding bound legs at the extremities of the cargo. Finally, the diffusivities in figure 2b are non-monotonic functions of the number of cargo legs. This is the result of competition between a total rate of events which increases with increasing N (leading to a higher rate of change in cargo position), and an average displacement per event which scales as ∼1/(n ± 1)² (see equations (2.1) and (2.2)).

In order to explore the effect of spatially varying cargo–substrate interactions, let us now consider a binding rate distribution of the form

for which we define the characteristic length scale l_c = 2(ln(39))^1/6 ≃ 2.5 as the distance from the origin to the position where κ₁(x) = 1.1κ₂. This distribution has a central region with higher binding rate than its surroundings (figure 3). As we can see in figure 3a, the result of a binding rate distribution given by equation (2.8) and a position-independent unbinding rate is that cargo exhibit an effective velocity profile that acts to accumulate the cargo near the centre of the high-binding-rate region. This is also true when ${\overline{k}}_3\mathrm{\Delta }x\ne 0$, in which case the stable fixed point for bound cargo motion is shifted in the direction of Δx, as predicted by equation (2.7) (see electronic supplementary material, figure S6). Figure 3 shows that the effective velocity and diffusivity calculated from equations (2.1), (2.2), (2.4) and (2.5) compare very well with the results from full numerical simulations. The small discrepancies are the result of the mean-field analytical approach using the steady-state probability of cargo having n legs bound, P_n(x), which is only correct for long-lived cargo. We find similarly good agreement between bound leg distributions from simulations and their corresponding heuristic fits (see electronic supplementary material, §6 and figure S4a).

• Download figure
• Open in new tab
• Download PowerPoint

Direct numerical simulations of cargo that interact with a substrate via the binding rate profile in equation (2.8) and a uniform unbinding rate (see figure 3) result in a probability distribution function P(x, t) that exhibits a strong central peak at sufficiently long times (figure 4a). The tracks of individual cargo clearly show that this is the result of the preferential binding of cargo legs within the central high-binding-rate region (figure 4b). Figure 4a also shows the probability distributions obtained by feeding the functions v_eff(x) and D_eff(x) (equations (2.4) and (2.5) and figure 3) into either the continuum model in equation (2.3), or stochastic molecular dynamics simulations of the Langevin dynamics in equation (2.6). These effective descriptions of cargo dynamics lead to probability distributions that agree remarkably well with the full model, while providing a simpler and more intuitive understanding of the behaviour of cargo. This is also true for the case where ${\overline{k}}_3\mathrm{\Delta }x\ne 0$, for which the distribution P(x, t) is skewed in the direction of Δx (electronic supplementary material, figure S6).

• Download figure
• Open in new tab
• Download PowerPoint

2.3. Behaviour of the effective velocity

When the cargo is sufficiently small compared to the characteristic lengths of variation in κ_1,2(x) (i.e. L ≪ |κ_1,2′(x_a(t))/κ_1,2″(x_a(t))|), the explicit dependence of $S_{\lambda }^{(1)}(x)$ on the underlying rate distributions can be simplified by using a leading order approximation to κ_1,2(x). Assuming a shape factor S_H( · ), this results in

where $I(x,n)=\int (y-x)P_b(y|x,n)\hspace{0.17em}\mathrm{d}y$ is the average difference between the unbinding position of a cargo leg and the cargo centre. We note that the distribution P_b(y|x, n) is symmetric only for uniform binding/unbinding rates, in which case $S_{\lambda }^{(1)}={\overline{k}}_3\mathrm{\Delta }x$. For non-uniform binding/unbinding rates, the first term in equation (2.9) predicts that the component of the effective velocity due to binding events is proportional to the local gradient in the binding rate distribution, whereas the second term states that the component due to unbinding events is dominated by variations in the bound leg distribution. For the cargo–substrate interaction rates defined in equation (2.8), it can be shown that the first term dominates over the second (see electronic supplementary material, figures S5 and S6c). In this case, equation (2.9) predicts, for ${\overline{k}}_3\mathrm{\Delta }x=0$, that $S_{\lambda }^{(1)}(x)\propto L^2{\kappa }_1^{\prime }(x)$, as well as that $S_{\lambda }^{(1)}(x)$ increases monotonically as a function of N until it plateaus at the value $\underset{N\to \mathrm{\infty }}{lim}(S_{\lambda }^{(1)}(x))=L^2\hspace{0.17em}{\kappa }_1^{\prime }(x)/3$ (see electronic supplementary material, §4). Equation (2.9) also predicts that cargo will accumulate around the maximum in the binding rate distribution. However, when ${\overline{k}}_3\mathrm{\Delta }x\ne 0$ the stable fixed point around which cargo accumulates shifts from the maximum of the binding rate towards one of the peaks in κ₁′(x), depending on the sign of Δx.

It is also possible for cargo to overhang the edges of a non-periodic substrate, in which case they will effectively experience step-like changes in their binding and unbinding rates (κ_1,2(x) = 0 beyond the edge) and L will no longer be small compared to the characteristic lengths of variation in κ_1,2(x) near the boundary. The dynamics of these cargo will still be described by equations (2.1), (2.2), (2.4) and (2.5), but no longer by equation (2.9). For binding/unbinding rates κ_1,2 which are uniform within the bounds of the substrate, the formalism in equations (2.1), (2.2) and (2.4) can be used to show that there will be a contribution to the effective velocity of the cargo of the form

where the choice of ± depends on whether we are considering the right (+) or left (−) edge of the substrate, resulting in $J_{+}(y|x,n,N)={\int }^{x_{\mathrm{edge}}}(y-x)P_b(y|x,n)\hspace{0.17em}\mathrm{d}y$ and $J_{-}(y|x,n,N)={\int }_{x_{\mathrm{edge}}}(y-x)P_b(y|x,n)\hspace{0.17em}\mathrm{d}y$. In this case, the bound leg distributions P_b(y|x, n) will be skewed in the direction of the substrate, meaning that the contribution to the effective velocity in equation (2.10) due to unbinding events will always act away from the substrate. By contrast, the binding event contribution to the effective velocity in equation (2.10) always acts towards the substrate. Since the effective velocity due to binding events is expected to dominate over the component resulting from unbinding events (see electronic supplementary material, figures S5 and S6c), cargo will exhibit an effective velocity that tries to maximize their overlap with the substrate. As a result, cargo can exhibit a stable fixed point in their motion near the edge of the substrate when ${\overline{k}}_3\mathrm{\Delta }x\ne 0$. This effect has been observed experimentally for actin filaments that interact with microtubules [5].

2.4. Using experimentally derived input parameters

Previous studies have shown that beads coated in EB binding domains can track the growing ends of microtubules in the presence of EBs [5,27], but a model that can describe these dynamics has not yet been developed. In order to test whether the model presented in this work can reproduce the dynamics of multivalent cargo in biologically relevant conditions, we first converted the parameters describing EB–microtubule interactions obtained in previous studies [23,27,36–41] into binding/unbinding rate distributions required as inputs for our stochastic cargo binding simulations (see electronic supplementary material, §7 and table S1). Here, we have subsumed the EB layer into an effective distribution of binding/unbinding rates for the multivalent cargo. Considering EBs explicitly does not fundamentally alter the results. Using these rate profiles, our model predicts that stable fixed points will be exhibited in the motion of (N ≥ 6)-legged cargo of sufficient size ($S_{\lambda }^{(1)}(x^{\ast })=0$), corresponding to them being able to track growing microtubule ends (figure 5). The component of each cargo’s effective velocity $S_{\lambda }^{(1)}(x)$ due to binding and unbinding events increases monotonically as a function of N until it plateaus, as predicted analytically in §2.3, and increases approximately quadratically as a function of L for comparatively small cargo, as predicted by equation (2.9) (see electronic supplementary material, figure S8).

• Download figure
• Open in new tab
• Download PowerPoint

2.5. Extending the cargo binding model to two dimensions

The cargo binding model defined in equations (2.1)–(2.5) can be extended to two dimensions to describe more complex cargo–substrate interaction networks. The PDF P(x, t) describing the probability of finding cargo at position x = (x, y) at time t can be derived using the same method as in the one-dimensional case (see electronic supplementary material, §8), and the resulting Fokker–Planck equation reads

andHere $S_{\lambda }^{(i,\cdot )}(\mathbf{x})$ and ${\lambda }_m^{(i,\cdot )}(\mathbf{x},n)$ are defined in the same way as their one-dimensional counterparts but with spatial integrals running over two dimensions, and $S_{\lambda }^{(i,xy)}(\mathbf{x})$ and ${\lambda }_m^{(i,xy)}(\mathbf{x},n)$ are functions of the two-dimensional cross moment $\int \mathrm{d}x\int \mathrm{d}y({[{\mathrm{\Delta }}_m^2(x-x_a(t))(y-y_a(t))]}^i{k}_m(\mathbf{x}|{\mathbf{x}}_a(t),n))$.

In a two-dimensional system, cargo legs exhibit the baseline interaction rates κ_m(x), which must be corrected using a two-dimensional shape factor S(x|x_a(t), n) and bound leg distribution P_b(x|x_a(t), n). Once again, the simple shape factor S_H(x|x_a(t)) can be introduced that limits the possible binding positions of cargo legs to within the range ∣∣x − x_a∣∣ ≤ L. Then, working through the same procedure followed for the one-dimensional case, we arrive at a complete set of equations for the two-dimensional model (see electronic supplementary material, §8). These equations show that two-dimensional cargo will move diffusively for position-independent binding/unbinding rates (see electronic supplementary material, figure S12a), while for spatially dependent rates they will exhibit an effective velocity towards regions of increasing binding rate and decreasing unbinding rate, as in the one-dimensional case. This can be seen, for example, in figure 6, which shows the effective velocity of cargo on a two-dimensional substrate with ${\overline{k}}_3=0$, a uniform unbinding rate κ₂, and a non-uniform binding rate defined by

which is the two-dimensional equivalent of equation (2.8). Similarly to what was observed in the one-dimensional case (figure 3a), we see that this binding rate distribution results in an effective velocity for bound cargo that acts purely towards the central region of increased binding rate.

• Download figure
• Open in new tab
• Download PowerPoint

2.6. Summary of key points

The model introduced here aims to provide a minimal conceptual framework to understand the noisy dynamics of a cargo interacting with a substrate through the binding and unbinding of its multiple interaction sites, here referred to as legs. We sought an approach that could connect these individual interaction events to the collective motion of the entire cargo. Cargo–substrate interactions begin with the binding of a single leg onto the substrate and end when the last remaining leg unbinds (figure 1). During this time, the position of the cargo is dictated by the average position of its bound legs (figure 1). For homogeneous substrates, this simple dynamics gives rise to an effective diffusion of the cargo with distinct short- and long-timescale diffusivities (figures 2b and 7a). The difference results from the time required for a cargo to equilibrate the initial distribution of its bound legs to the long-timescale one (see §2.2). We propose that this could explain previous reports of decaying cargo diffusivity at long timescales [5].

• Download figure
• Open in new tab
• Download PowerPoint

Substrates with spatially varying binding/unbinding rates, instead, will give rise to an effective diffusivity D_eff(x) and an effective velocity v_eff(x) dependent on position. These can be connected to the binding/unbinding rate distributions through the analytical expressions in equations (2.4) and (2.5) (figure 7b). The drift is caused by differences in the binding/unbinding rates experienced by different parts of a cargo, leading it to accumulate towards regions of higher binding rates and/or lower unbinding rates (figure 4a). Differences in binding and unbinding rates can alternatively be understood as differences in a putative free energy landscape, such that higher binding rates correspond to lower free energies and a stronger ‘attraction’. For a fixed substrate, the system then evolves towards a final steady-state distribution akin to that of a free energy minimum (see for example figure 4a), in a way that is conceptually similar to a Brownian particle drifting down a potential well. However, for substrates that grow continuously such as a polymerizing microtubule, the system is kept far from equilibrium by the constant addition of new binding sites and the maturation of old ones (see electronic supplementary material, §7). In this case, if the spatial gradient of binding/unbinding rates is sufficient given the size of the cargo and the number of its legs (see §2.3), the cargo will be able to ‘surf’ the wave of newly added binding sites and follow the growth of the substrate (see figure 7c and §2.4). In the Brownian particle analogy, this is equivalent to the potential well exhibiting a constant drift, which causes the equilibrium position of the Brownian particle to be displaced from the centre of the well.

4.1. Simulation methods

Stochastic cargo binding simulations were implemented in Matlab using the Gillespie algorithm [46,47] to probe the system state in continuous time (see electronic supplementary material, §9). Molecular dynamics simulations were also implemented in Matlab, but instead updated cargo positions according to the Langevin equation defined in equation (2.6) using a forwards Euler scheme. Wiener process displacements were calculated using inverse transform sampling and binding dynamics were introduced by randomizing the position of the cargo within the periodic domain with rate k_ran(x) = 1/t_dwell(x) (calculated by substituting the local binding and unbinding rates into a previously published average dwell time formula [29,30]).

4.2. Numerical methods

Numerical solutions to the Fokker–Planck equation defined in equation (2.3) were obtained using the built-in Matlab function pdepe().