## Abstract

This paper presents an analytical solution of one-dimensional transient molecular-motor-assisted transport equations that describe transport of different organelles (such as transport vesicles loaded with a cargo of specific proteins) in a neuron's axon or dendrite. Large intracellular organelles are transported in the cytoplasm by a combined action of diffusion and motor-driven transport. In an axon, organelles are transported away from the neuron's body towards the axon's terminal by kinesin-family molecular motors running on tracks composed of microtubules (MTs); old and used components are carried back towards the neuron's body by dynein-family molecular motors. Using the method of separation of variables, a generalized Fourier series solution for this problem is obtained. The solution uses three different orthogonal sets of eigenfunctions to represent the concentration of free organelles transported by diffusion, MT-bound organelles transported away from the neuron's body, and MT-bound organelles transported towards the neuron's body. Binding/detachment kinetic processes between the organelles and the MT are specified by first-order rate constants; these lead to coupling between the three organelle concentrations.

## 1. Introduction

Diffusion is not a sufficiently fast mechanism for transporting large intracellular particles (organelles), such as large protein particles or intracellular vesicles, carrying different types of cargo between intracellular compartments (for example, between the endoplasmic reticulum, where most protein particles are synthesized, and the Golgi apparatus, where they are sorted and dispatched (Alberts *et al*. 2002)). This is because according to Einstein's relation that determines the diffusivity of small particles due to the Brownian motion, the diffusivity is inversely proportional to the particles' radius, which means that larger particles have smaller diffusivity. To overcome the limitations of diffusion transport, cells rely on directional transport of organelles along the cytoskeleton that is composed of different types of protein filaments such as microtubules (MTs); this transport is powered by molecular motors. Depending on the type of molecular motor attached to the intracellular cargo, it is transported towards either the plus or minus end of a MT. Recently, Smith & Simmons (2001) have suggested equations governing molecular-motor-assisted transport of intracellular organelles. The model suggested in Smith & Simmons (2001) describes the motion of intracellular particles under the combined action of diffusion and motor-driven transport. The latter occurs on a network of intracellular filaments (MTs or actin cables) that provide a stable framework along which molecular motors can generate the force necessary to move the particle along the MT. According to this model, an organelle attaches itself to a molecular motor (a specialized protein that, as a result of a chemical process, usually ATP hydrolysis, undergoes a conformational change ‘walking’ along intracellular filaments, such as MTs). The organelle either diffuses freely in the cytosol or moves on a filament at a motor velocity, and can bind to or detach from a filament. Depending on the type of a molecular motor (or several motors) attached to the particle, the motion along the MT can occur in either direction.

Dinh *et al*. (2005) presented a numerical solution of the Smith–Simmons equations to describe intracellular trafficking of adenoviruses between the cell membrane and the cell nucleus. Other aspects of the intracellular transport of cell organelles and vesicles along MTs are considered in Ajdari (1995), Lipowsky *et al*. (2001), Nieuwenhuizen *et al*. (2002), Welte (2004), Friedman & Craciun (2005), Pangarkar *et al*. (2005), Dinh *et al*. (2006), Leopold & Pfister (2006), Kuznetsov (in press *a*,*b*) and Kuznetsov *et al*. (in press).

One-dimensional Smith–Simmons equations in Cartesian coordinates can be used to describe the transport of organelles in long specialized arms (processes) of neurons (figure 1). If the arm transmits electrical signals it is called an axon and if it receives electrical signals it is called a dendrite. Axons in a human body can be up to 1 m in length. All MT in an axon have the same polarity (their plus ends point towards the axon terminal); they do not stretch the entire length of the axon and the continuous path along the axon is composed of short overlapping segments of parallel MTs. Transport vesicles loaded with specific proteins are carried away from the neuron body towards the synapse (the axon terminal) by kinesin-family molecular motors (this family of molecular motors is responsible for the transport on MTs towards their plus ends). Used and old components are carried from the axon terminal towards the body of the neuron by dynein-family molecular motors. (This family of molecular motors is responsible for the transport on MTs towards their minus ends.) In dendrites the MT polarities are mixed; some of them point their plus ends towards the dendrite tip and some towards the neurons' body. Therefore, in a dendrite, depending on the polarity of a particular MT, transport in a certain direction (to the neuron's body or away from it) can be carried out by either kinesin or dynein molecular motors (Alberts *et al*. 2002; Pollard & Earnshaw 2002). Irregularities in axonal transport may lead to traffic jams resulting in axonal swellings that contain abnormal amounts of molecular-motor proteins and various organelles; these traffic jams have been linked to neurodegenerative diseases such as Alzheimer's disease (Hurd & Saxton 1996; Goldstein 2001; Stokin *et al*. 2005).

Kuznetsov (2007) obtained an analytical solution of the Smith–Simmons equations through generalized Fourier series using a set of non-orthogonal eigenfunctions. In this paper, the governing equations suggested in Smith & Simmons (2001) are extended by accounting for the random component of motion of the intracellular organelles bound to MTs. The locomotion generated by molecular motors can be split into (average) deterministic and random components. The randomness in the locomotion is caused by different reasons; one is that the molecular motor works in a very noisy environment constantly experiencing thermally excited collisions with water molecules (Lipowsky *et al*. 2006). There are several approaches to modelling the stochastic component of motion induced by molecular motors (for example, this can be done by calculating the time evolution of probabilities of different motor states); one popular approach is based on the assumption that the random component can be approximated by a diffusive process. Using this assumption, molecular-motor-assisted transport equations are extended; an analytical solution of one-dimensional transient molecular-motor-assisted transport equations is obtained through generalized Fourier series using three orthogonal sets of eigenfunctions following the approach outlined in Cain & Meyer (2006).

## 2. Governing equations

The one-dimensional model suggested in Smith & Simmons (2001) describes the transport between two planar boundaries, and , where all the particle concentrations vary in the -direction only (figure 1). A fraction of space between the boundaries is occupied by MTs (which are assumed to be strictly parallel to the -direction); the remaining space allows for the diffusion of unbound (free) particles in the -direction. In an axon (figure 1), if a particle bound to a kinesin-family molecular motor attaches to a MT, the force generated by this molecular motor transports this particle with a positive velocity, , away from the neuron cell body towards the axon terminal. If a particle bound to a dynein-family molecular motor attaches to a MT, the force generated by this molecular motor transports this particle with a negative velocity, , away from the axon terminal towards the neuron cell body. Unbound particles may attach to MTs; corresponding first-order rate constants for binding to MTs for the particles moving in the positive and negative directions are and , respectively. The particles attached to MTs may detach from them; corresponding first-order rate constants for detachment from MTs for the particles moving in the positive and negative directions are and , respectively (owing to the thermal noise, an average molecular motor falls off from a MT after it has made approximately 100 steps (Lipowsky *et al*. 2006)). The governing equations describing particle transport in a neuron's arm (a dendrite or axon) presented in Smith & Simmons (2001) are extended by assuming that the random component of motion for the particles bound to MTs can be approximated by a diffusive process. The random component of the motor movement is typically small (Visscher *et al*. 1999); however, the randomness can be caused by obstacles in either the cytoplasm or the filament. Another source of randomness would be short backward movements or pauses either inherent to dynein or dynein-like motors (Ross *et al*. 2006), or the presence of motors on the cargo that pull in the opposite direction (Muller *et al*. 2008). The dimensional forms of governing equations are as follows:(2.1)(2.2)(2.3)where is the diffusivity of a free particle; are the diffusivities of MT-bound particles, modelling random aspects of molecular-motor-assisted transport; is the time; is the free-particle concentration; is the concentration of particles moving on MTs in the positive direction (away from the cell body); is the concentration of particles moving on MTs in the negative direction (towards the cell body); is the linear coordinate along the axon; is the velocity of a particle moving on a MT towards the cell body (in an axon this is the motor velocity generated by a dynein-family molecular motor), is negative; is the velocity of a particle moving on a MT away from the cell body (in an axon this is the motor velocity generated by a kinesin-family molecular motor), is positive; and are the first-order rate constants for binding to MTs for the particles that move in the positive and negative directions, respectively; and and are the first-order rate constants for detachment from MTs for the particles that move in the positive and negative directions, respectively. For , equations (2.1)–(2.3) collapse to equations suggested in Smith & Simmons (2001). There is an additional benefit of adding the diffusion terms to equations (2.2) and (2.3). Without these terms, equations (2.2) and (2.3) are wave equations whose numerical solutions, unless special numerical schemes are used, exhibit oscillatory behaviours. The addition of diffusion terms to equations (2.2) and (2.3) removes these unphysical oscillations from the numerical solution.

Since the model is formulated in terms of three concentrations, the concentration of free particles, , particles bound to MTs moving away from the neuron's body, , and particles bound to MTs moving towards the neuron's body, , the model presented here is not restricted to any particular type of molecular motors responsible for the transport on MTs, and equations (2.1)–(2.3) are applicable not only to organelle transport in axons, but also in dendrites.

Equations (2.1)–(2.3) are solved subject to the following boundary conditions:(2.4)(2.5)where and are fixed concentrations of particles at and , respectively, and *σ*_{0} and *σ*_{L} are the degrees of loading at and , respectively.

The initial condition is(2.6)

## 3. Dimensionless equations

According to the Buckingham Pi theorem (White 2008), the maximum reduction is equal to 2 (the number of dimensions describing the variables, length and time). Dimensionless variables and dimensionless groups (parameters) are introduced as follows:(3.1)(3.2)(3.3)The dimensionless governing equations are(3.4)(3.5)(3.6)The dimensionless boundary conditions are(3.7a–c)(3.8a–c)The dimensionless initial condition is(3.8d)

## 4. Solution technique

Because equations (3.4)–(3.6) are linear, the problem is split into a steady non-homogeneous problem and an unsteady homogeneous problem as follows:(4.1)where , and satisfy the following equations:(4.2)(4.3)(4.4)Equations (4.2)–(4.4) must be solved subject to the following boundary conditions:(4.5a–c)(4.6a–c)The exact solution of steady-state equations (4.2)–(4.4) with boundary conditions (4.5*a–c*) and (4.6*a–c*) is given in appendix A.

The transient components, , and , satisfy the following equations:(4.7)(4.8)(4.9)Equations (4.7)–(4.9) must be solved subject to the following boundary conditions:(4.10a–c)(4.11a–c)and the following initial condition:(4.12a–c)

The following linear operators are introduced:(4.13)(4.14)(4.15)Using definitions of *L*_{1}, *L*_{2}, and *L*_{3}, equations (4.7)–(4.9) are recast as(4.16)(4.17)(4.18)The eigenvalue problem associated with *L*_{1} and the boundary conditions for are(4.19)(4.20)(4.21)The solution of equations (4.19)–(4.21) gives the following orthogonal set of eigenfunctions:(4.22)The eigenvalue problem associated with *L*_{2} and the boundary conditions for is(4.23)(4.24)(4.25)The solution of equations (4.23)–(4.25) gives the following orthogonal set of eigenfunctions:(4.26)The eigenvalue problem associated with *L*_{3} and the boundary conditions for is(4.27)(4.28)(4.29)The solution of equations (4.27)–(4.29) gives the following orthogonal set of eigenfunctions:(4.30)It is assumed that the functions , and can be represented by the following generalized Fourier series:(4.31)(4.32)(4.33)where , and are the orthogonal eigenfunctions given by equations (4.22), (4.26) and (4.30), respectively. Generally speaking, sums in the series in equations (4.31)–(4.33) are infinite sums (*N*→∞), but for practical reasons their evaluation is approximated by a finite sum within the required limits of desired accuracy; the convergence of these series is demonstrated in the figures presented later on.

Substituting equations (4.31)–(4.33) into equation (4.16), the following is obtained:(4.34)Since eigenfunctions are orthogonal, multiplying equation (4.34) by and integrating results in(4.35)Equations for the functions *α*_{1}(*t*), *α*_{2}(*t*), *α*_{3}(*t*) and *α*_{4}(*t*) for *N*=4 are given in appendix B.

The initial conditions for functions *α*_{l}(*t*) [*l*=1, …, *N*] are obtained by projecting the difference between the initial and steady-state distributions on the subspace composed of first *N* orthogonal eigenfunctions, (*i*=1, …, *N*),(4.36)Multiplying equation (4.36) by and integrating results in(4.37)Equations for the initial values *α*_{1}(0), *α*_{2}(0), *α*_{3}(0) and *α*_{4}(0) for *N*=4 are given in appendix C.

Substituting equations (4.31)–(4.33) into equation (4.17), the following is obtained:(4.38)Since eigenfunctions are orthogonal, multiplying equation (4.38) by and integrating results in(4.39)The initial conditions for the functions *β*_{h}(*t*) [*h*=1, …, *N*] are found from(4.40)Multiplying equation (4.40) by and integrating results in(4.41)Substituting equations (4.31)–(4.33) into equation (4.18), the following is obtained:(4.42)Since eigenfunctions are orthogonal, multiplying equation (4.42) by and integrating results in(4.43)The initial conditions for the functions *γ*_{f}(*t*) [*f*=1, …, *N*] are found from(4.44)Multiplying equation (4.44) by and integrating results in(4.45)equations (4.35), (4.39) and (4.43) give 3*N* first-order ordinary differential equations for the functions *α*_{l}(*t*) [*l*=1, …, *N*], *β*_{h}(*t*) [*h*=1, …, *N*] and *γ*_{f}(*t*) [*f*=1, …, *N*], respectively, and initial conditions for these equations are given by (4.37), (4.41) and (4.45), respectively. Once these functions are determined, the solution for , and is given by equations (4.31)–(4.33).

## 5. Results and discussion

Computations are performed for the values of dimensionless parameters summarized in table 1. Figures 2–4 demonstrate a good convergence of the generalized Fourier series solution to the numerical solution (obtained with very small Δ*t* and Δ*x* to ensure high accuracy) for three dimensionless times, *t*=2,4 and 8. These figures demonstrate that the convergence is quite fast, and that the generalized Fourier series solution given by equations (4.31)–(4.33), truncated to the first six terms, is within a few per cent of the high-accuracy numerical solution.

Figure 5 shows distributions of the dimensionless concentrations of intracellular particles (*n*_{0}, *n*_{+} and *n*_{−}) at different times (*t*=0 (initial condition), 2,4,8 and steady state) obtained numerically by a high-accuracy numerical solution. The behaviour of the free-particle concentration, *n*_{0}, is quite interesting (figure 5*a*). The initial concentration of *n*_{0} is 0.1. Then, at *t*=2, the concentration of free particles falls significantly below the steady-state concentration. As time increases, *n*_{0} increases and at *t*=8 it almost reaches the steady-state concentration. This effect of overcompensation (which is similar to overstability) is apparently related to the coupling of *n*_{0} with two other concentrations, *n*_{+} and *n*_{−}, which initially have uniform zero concentrations. Owing to that at the beginning of the process free particles intensively bind to MTs, which causes a fast decrease of free-particle concentration. As this happens, the concentration of free particles falls below the steady-state concentration and then it grows back to the steady-state concentration as time increases.

It is also interesting that for *t*=8, both *n*_{0} and *n*_{+} have almost reached their steady-state concentrations (figure 5*a*,*b*), but *n*_{−} is still quite far from steady state (figure 5*c*). The total concentration *n*_{0}+*n*_{+}+*n*_{−} is the quantity accessible to experiments; its distribution at different times is shown in figure 5*d*. This figure shows that the total concentration of organelles increases as time increases.

## 6. Conclusions

The exact solution of the one-dimensional steady molecular-motor-assisted transport equations, that describe transport of different organelles in a neuron's axon or dendrite, is obtained. The transient solution of this problem is obtained through generalized Fourier series using three different orthogonal sets of eigenfunctions. The series solution is demonstrated to converge fast to a high-accuracy numerical solution. The analysis of the transient behaviour of the free-particle concentration reveals that it may drop to a low concentration, passing the steady-state concentration, and then grow back to the steady-state concentration. This effect of overcompensation is explained by zero initial concentrations of particles bound to MTs used in the presented computation, which causes a fast initial drop of the free-particle concentration as free particles bind to MTs.

## Acknowledgments

The authors gratefully acknowledge NATO Collaborative Linkage Grant (CBP.NUKR.CLG 981714). The authors are indebted to the reviewers whose critical comments helped improve this manuscript.

## Footnotes

- Received March 25, 2008.
- Accepted May 16, 2008.

- © 2008 The Royal Society