## Abstract

We present a computer simulation study, via lattice Boltzmann simulations, of a microscopic model for cytoplasmic streaming in algal cells such as those of *Chara corallina*. We modelled myosin motors tracking along actin lanes as spheres undergoing directed motion along fixed lines. The sphere dimension takes into account the fact that motors drag vesicles or other organelles, and, unlike previous work, we model the boundary close to which the motors move as walls with a finite slip layer. By using realistic parameter values for actin lane and myosin density, as well as for endoplasmic and vacuole viscosity and the slip layer close to the wall, we find that this simplified view, which does not rely on any coupling between motors, cytoplasm and vacuole other than that provided by viscous Stokes flow, is enough to account for the observed magnitude of streaming velocities in intracellular fluid in living plant cells.

## 1. Introduction

Within eukaryotic cells, intracellular flows are often unremarkable: typical Reynolds numbers are low, and cells are tiny enough to quench most fluid motion within them. An important exception is provided by cytoplasmic streaming, which is the name given to the approximately 80–100 µm s^{−1} directed flow of cytosol and organelles around large plant and fungal cells [1–3].^{1} Unlike typical mammalian cells, which are tens of micrometres in size, the internode stalks of the *Chara corallina* algae, for instance, are single cells that are about a millimetre in diameter and several centimetres long. Giant single cells such as these need to overcome a highly non-trivial transport problem in order to deliver nutrients and other chemicals throughout their interior. Given such length scales, thermal diffusion does not provide a viable option to move things around quickly enough: for instance, given typical intracellular diffusion coefficients, it would take months for even highly mobile ions to traverse the cell length. Cytoplasmic streaming gives a very efficient alternative, as it sets the whole fluid within the cell in motion, thereby readily advecting proteins and nutrients which happen to be there—in other words while the Reynolds number may be infinitesimally small, the Peclet number is instead rather large [5,6].

Understanding the basic biophysical mechanism underpinning cytoplasmic streaming has been the subject of an ongoing debate at the interface between fluid dynamics and cell biology. To introduce the competing scenarios proposed in the literature, it is useful to recall the geometry of algal cells such as *Chara corallina*. The basic components are the ‘internodes’ or internode stalks, which are extra-long single cells as mentioned earlier. Each of these internodes has an approximately cylindrical symmetry, and may be thought of as consisting of two concentric layers: the external one is a micrometre-thick cellular layer, filled with cytosol (endoplasm) and separated from the large interior, known as the vacuole, by a fluid membrane. The vacuole is made up of an aqueous solvent, and contains no internal structure. The inner surface of the external wall of the internode cell is instead covered by actin bundles, which act as conveyor belts along which myosin motors walk. The actin bundles, or lanes, follow helical paths on the cylindrical walls of *Chara*, and they are organized into alternating bands having opposite directions of velocity. The bands are separated by a so-called indifferent zone. The hypothesized mechanism for cytoplasmic streaming is that the directed motion of myosin motors on the actin bundle tracks can somehow entrain the cytoplasmic fluid, and that this in turn sets the vacuole into motion by transferring momentum through the fluid membrane separating it from the cytoplasm. According to this view, the myosin traffic provides an effective shear boundary condition on a Stokes equation for the velocity inside the vacuole. The fact that ‘the tracks along which motors move are helicoidal’ leads to an interesting fluid dynamics problem. Its solution, presented in van de Meent *et al.* [5], includes a secondary flow along the planes perpendicular to the cylindrical axis, which helps nutrient and chemical mixing along the flow gradient direction as well—mixing is a notoriously hard requirement to fulfil at zero Reynolds number. More recently, direct nuclear magnetic resonance velocimetry experiments of flow inside the vacuole of single internodal cells have quantitatively validated this analytical solution [7].

Earlier experiments on cytoplasmic streaming based on laser light scattering experiments on *Nitella* [8,9] and *Chara* [10] cells showed that the fluid velocity distribution is quite narrow. While Sattelle & Buchan [10] speak of plug-like motion for both *Chara* and *Nitella*, Mustacich & Ware [9] emphasize the fact that the peak in the spectrum of scattered light is not sharp enough for the streaming flow of a plug but indicates that most particle velocities are within 10 per cent of the most likely velocity. Experiments support that it is unlikely that any particles move at substantially larger velocities than those carried by cytoplasmic streaming—in other words, the fluid velocity has to be very close to that of the vesicles dragged along by the motors and driving the flow.

The idea that myosin motion can provide an effectively smooth boundary condition for the fluid dynamics inside the vacuole by necessity requires a strong hydrodynamic coupling between actin–myosin motion on the cell surface and fluid flow in the vacuole. Finding exactly what structure can provide this coupling is however a non-trivial task. In an interesting contribution, Nothnagel & Webb [11] analysed the hydrodynamic feasibility of three models for momentum transfer from myosin to endoplasm and vacuole. Their calculation shows that individual myosin molecules running on the actin tracks are by themselves ineffective in setting the cytosol or the vacuole into motion. A second model considered attachment of myosin to organelles and vesicles in the endoplasm: while this much improves the viscous coupling, the calculations in the study of Nothnagel & Webb [11] still suggested that if this mechanism was at the origin of cytoplasmic streaming it would require very packed traffic of motors along the actin lanes, which is unlikely in the real system. The final model, favoured by the semianalytical estimates in the work of Nothnagel & Webb [11], envisages the existence of an elastic network, or gel, incorporating the moving motors and extending into the endoplasm: in this framework, the movement of the motors pulls the network forward in a plug-like fashion and the coupling to the vacuole is readily achieved.

Subsequent experiments (video-enhanced light microscopy and fast-freezing electron microscopy) further suggested that the myosin motors may actually be attached to the endoplasmic reticulum, which then performs a sliding motion over the actin cables [12]. In some later papers, this view has been adopted [13], while others are less definite and speak only of cargoes in general [14] or state that the organelles associated with myosins have not been identified [2].

Here, we revisit the view that the myosin motors need not be directly associated with a network structure but that the high viscosity of the cytoplasm together with a thin slip layer are sufficient to make it move with a very flat velocity gradient at roughly the same velocity as the active particles. We perform lattice Boltzmann (LB) simulations of micrometre-sized spherical vesicles moving close to a planar wall (dragged by the motors) and show that, depending on the slip allowed at the wall, this simpler picture can explain the uniform streaming profile also for realistic densities of spherical particles. The effect of the endoplasmic reticulum and all other cell contents is incorporated into the high viscosity of the cytoplasm, but our results suggest that it is not necessary to invoke a physical tethering of the motors to any network in order to explain the known phenomenology of cytoplasm streaming.

## 2. Model and methods

### 2.1. Lattice Boltzmann simulations of colloidal spheres in a fluid

Rather than address the cylindrical geometry suggested by the *Chara* system, we address here a simpler, planar model (figure 1). This allows us to address the issues of principle outlined earlier while avoiding the numerical difficulties of a curved geometry. Moreover, as detailed later, the physics of interest to us concerns the cytosol and only a thin region of the vacuole containing the endoplasm–vacuole interface: this region is almost planar. We consider a slab of material made up of two layers with different viscosities, *η*_{I} and *η*_{II} > *η*_{I}, corresponding to, respectively, the vacuole and the endoplasm, modelled here as a Newtonian fluid with high viscosity, as in the semianalytical treatment in Nothnagel & Webb [11]. In each of the layers, the fluid obeys the Navier–Stokes equation,
2.1where *ρ* denotes fluid density, **u** is fluid velocity, *p* is pressure, *η*_{i} is the fluid viscosity and *i* = *I*,*II* labels the layer under consideration. In our formalism, Greek letters denote Cartesian components and summation over repeated indices is implied. In order to solve equation (2.1), we employ our Ludwig code, which is based on an LB algorithm (see Cates *et al.* [15] for details). Briefly, the LB method is based on a set of mesoscopic velocity distribution functions, *f*_{i}(**x**,*t*), which depend on some lattice velocity vectors **e**_{i}. In our three-dimensional case, we consider 19 velocity vectors (see appendix A); this is known as a D3Q19 lattice. The density and velocity fields of the fluid may be recovered as moments of such distribution functions, as follows:
2.2

In our LB algorithm, the distribution functions evolve according to the following discretized, or lattice, Boltzmann equation:
2.3where Δ*t* is the time step (normally unity in LB), *τ* is a relaxation parameter related to the fluid viscosity *η* via *η* = (*ρ*/3)(*τ* − Δ*t*/2), and *f*_{i}^{eq} is a set of equilibrium distribution functions. Note that *τ* differs in the vacuole and cytoplasm as the viscosities differ. The macroscopic equations of motion obeyed by *ρ* and *u*_{α} may be found by taking moments of equation (2.3) and expanding for small Δ*t*—a formal procedure known as the Chapman–Enskog expansion. The LB method defines some dynamical rules for the distribution functions that correspond to a continuity equation for the fluid density and equation (2.1) for its velocity field. It can be shown that this is the case, provided the following constraints on the moments of the equilibrium distribution are satisfied:
2.4

At this point, it should be noted that any set of equilibrium distribution functions satisfying these constraints leads to a viable LB algorithm to solve the fluid dynamics problem we are interested in. A popular choice in the literature is to expand the equilibrium distribution functions as a power series in the velocity vectors, with the coefficients being determined, in general in a non-unique way, through the constraints. Note that as we want to model a two-layered fluid, made up by a more viscous endoplasm and a less viscous vacuole, we considered a generalization of the LB equation given earlier, to allow for a spatially dependent value of the relaxation constant *τ*. We found that this simple procedure describes our two-layered system well, but we anticipate that this may not hold in more complicated systems where the interface is not planar.

The LB method, in the formulation which we adopt, has been developed to allow consistent coupling between the fluid dynamics of a Newtonian fluid and the dynamics of spherical colloidal particles embedded in it. The coupling is through the well-known method of bounce-back on links, which accurately enforces no-slip boundary conditions for the velocity field on the surface of each particle. Instead of repeating the discussion of these colloidal boundary conditions here, we refer the interested reader to Cates *et al.* [15] and Nguyen & Ladd [16] for details. In our case, the solid particles represent the motors with their associated cargoes (vesicles and organelles). They are positioned at a fixed distance from the bottom wall, and are subjected to a strong harmonic potential that virtually confines their motion to fixed one-dimensional lines, which represent the actin bundles.

A noteworthy advantage of the LB methodology we use is that, as it fully handles the boundary conditions for the velocity fields both on the particle surface and on the walls, both near- and far-field effects are correctly included. With respect to, for example, Stokesian dynamics, which would be another perfectly valid choice for the problem at hand, LB requires one to introduce both an inertial term and a finite compressibility: both these terms need to be kept small so that they do not effectively contribute to the physics of our system. On the other hand, LB allows one to reach much larger systems than Stokesian dynamics, at least for the non-uniform geometry addressed here. Finally, we note that an analytical approximation based on modelling the vesicles as Stokeslets near a wall would not incorporate near-field interactions that are important in our simulations, as particles may come into close contact.

### 2.2. Boundary conditions for slip and no-slip walls

In this section, we discuss in some detail how to implement slip boundary conditions at the bottom wall (see geometry in figure 1), which is important for our work. Here, we outline the main equations that we use in our algorithm, whereas a full derivation is given in appendix A.

In order to model slip at a solid wall in LB simulations, a partial slip parameter *p* can be introduced that intuitively corresponds to the ‘fraction of slip’ at the wall. More accurately, *p* is the fraction of each of the distribution functions that is ‘bounced forward’ at the wall, whereas 1 − *p* is the fraction that is ‘bounced back’ (see appendix A). The limit of *p* → 0 corresponds to the commonly employed no-slip boundary condition. The amount of wall slip may also be related to the Knudsen number of the fluid [17]. In our scheme, we use it as a parameter to be determined, in practice, from the knowledge of the slip length, *u*(0)/*u*′(0), which we assume comes from either experiment or theory.

Our choice of a slip wall for the bottom boundary is motivated by the need to model the fluid dynamics in the endoplasm close to the dragged cargo. We assume that below the cargo lies a thin aqueous low-viscosity layer, as the viscosity in the endoplasm derives either from the presence of a labile reticular network with mesh size smaller than the distance between sphere and wall or from macromolecular crowding in the endoplasm. Either way, the network or crowding agents, which are coarse-grained into a high-viscosity fluid for the bulk endoplasm, are depleted in the vicinity of the wall, and this accounts for the thin layer of lower viscosity. If the behaviour of the thin layer itself is not of interest, all it does is to introduce effective slip for the main system, which can be simulated using the slip parameter *p*.

To make progress, we start by writing down the velocity profile of the fluid as a power series in the distance from the bottom wall *z*,

The slip length *u*(0)/*u*′(0) can be determined analytically for *ɛ* ≪ *R*, where *R* is the total system size and *ɛ* is the thickness of the low-viscosity layer. By enforcing no slip at the bottom wall, together with continuity of the stress and velocity fields at *z* = *ɛ*, and subsequently taking the limit *ɛ* → 0 (see appendix A), we obtain the form for the slip length,
2.5to first order in *ɛ*, where *η*_{I} is the (lower) viscosity of the thin layer close to the wall and *η*_{II} is the (higher) viscosity in the bulk of the system. Here, we assume that the viscosity of the thin layer is the same as that of the vacuole as both are basically aqueous fluids. This is however only a conceptual simplification and not necessary for the derivation of the boundary condition.

The slip length *u*(0)/*u*′(0) and slip parameter *p* may be shown to be related by the following formula (derived in appendix A):
2.6

The last two equations also lead to the following explicit formula for *p*:
2.7

As demonstrated in appendix A, numerical simulations of simple test cases performed with a given *p* indeed lead to a flow profile that agrees with our analytical calculation. In the formulae given earlier, all quantities are given in terms of lattice, or simulation, units. To convert these to physical units, *ɛ* has to be multiplied by the lattice spacing Δ*x*, while *η* should be multiplied by the fluid density *ρ* and by Δ*x*^{2}, and divided by the time step Δ*t*. In lattice, or simulation, units, *ρ* = Δ*x* = Δ*t* = 1.

In physical units, we choose Δ*x* = 0.4 µm, Δ*t* = 0.08 ms and *ρ* ∼ 2.5 × 10^{7} kg m^{−3}. With these values, our vesicle size and vacuole viscosity map to 0.5 µm and 1cP, respectively, as in the experiments. Furthermore, a velocity of 0.01 in LU, typically used in our simulations, maps to a velocity of 50 µm s^{−1}. However, our chosen density is much higher than that of water.^{2}

### 2.3. Choice of parameters

Simulations were run for colloidal particles of radius *a*_{0} = 1.25 LU, whereas the size of the cytoplasm was varied within a physically meaningful regime. In keeping with the parameters of Nothnagel & Webb [11], vesicles are assumed to have a radius of 0.5 µm, meaning that 1 LU = 0.4 µm (see above). The spacing between actin cables (or ‘lanes’) was fixed to *d*_{l–l} = 2 µm = 5 LU. For computational reasons, it is desirable for the distance of particles from the wall to be such that there is at least one fluid node between the wall and the particle surface (otherwise, lubrication forces need to be considered). Particle centres are therefore placed at *z*_{p} = 2.35 LU (and they do not extend below *z* = 1.1 LU). The distance between the wall (mid-grid at *z*_{w} = 0.5) and the particle surface is thus *d*_{w–p} = 0.6 LU. This agrees well with the actual distance in the biological system: diameters of actin cables are between 100 and 200 nm [19], and the size of *Chara* myosin is estimated to be about 100 nm [14], much larger than typical animal myosin motors. This results in about *d*_{w–p} = 0.2 − 0.3 µm = 0.5 − 0.75 LU. Initial spacing between particles (centre to centre) varies with density: a separation *d*_{p–p} = 12.5 LU for instance corresponds to a line density of 0.2 (in units where 1 corresponds to a solid line of particles within a lane).

We fixed the vacuole viscosity to *η*_{I} = 0.02 and varied *η*_{II} in line with existing estimates—see Nothnagel & Webb [11], who assumed a ratio *η*_{II}/*η*_{I} of between 6 and 250. In simulations, the vacuole has thickness 75 LU = 30 µm. While this is much too small as the vacuole thickness can be approximately 500 µm, the location of the upper boundary has virtually no effect on the velocities on the endoplasm–vacuole interface, which we are ultimately interested in. We also note that moving the upper vacuole boundary further away would, if anything, only increase the cytoplasm's velocity with respect to particle velocities.

Finally, we note that here we assume that the viscosity of the thin layer in the cytoplasm that cannot be reached by the crowding agents is the same as that of the vacuole. This results in a slip layer of thickness equal to 220 nm, when *p* = 0.9. The thickness of the low-viscosity layer is not known but is restricted by the particle-to-wall distance (*d*_{w–p}) being equal to 0.6 LU, because the particles have to move within the high-viscosity region to provide effective drag on the main endoplasm layer.

## 3. Results

Figure 1 shows the geometry of our simulations, which we recall consists in a two-layer fluid, modelling the endoplasm–vacuole system, with a slip length boundary condition (equation (2.6)). A two-dimensional carpet of molecular motors with their attached cargoes is arranged along regularly spaced actin cables and each of the motors is further subject to a constant external forcing that drives it along the lane.

### 3.1. Allowing for a finite wall slip can reproduce the phenomenology of cytoplasmic streaming

Figure 2 shows the mean fluid velocity *u*(*z*) as a function of *z*, averaged over *x* and *y* and normalized by the average velocities of the motor carpet—the latter are recomputed for each simulation as they may vary. Figure 2*a* considers no-slip boundary conditions, *p* = 0. Our direct simulations of the three-dimensional flow within the endoplasm suggest that the velocities within the fluid sharply drop to about 50 per cent of the motor–cargo velocities close to the wall. Interestingly, this is very much in line with the semianalytical calculations presented in the study of Nothnagel & Webb [11] and relying on more assumptions than our LB simulations. As observed in the work by Nothnagel & Webb [11], this is not compatible with experimental measurements that suggest that the velocity distribution within the endoplasm is within 10 per cent of the motor velocities on the actin cables. The no-slip simulations also show a small peak in *u*(*z*), just above the particles. If the size of the particles is increased, such a peak disappears, so this is probably a discretization artefact.

Figure 2*b* instead shows what happens if we allow for wall slip. As mentioned in §2, in the simplest approximation, the slip may be characterized by a single length scale, henceforth denoted as the slip layer thickness. This is equal to 220 nm in figure 2*b*, which is a reasonable value given the size of the cargoes, and corresponds to a slip parameter *p* = 0.9 in our LB simulations. The cytoplasmic flow is now markedly different and much faster with respect to the no-slip case, and a tracer would move at over 90 per cent of the mean motor speed, in agreement with the observations.

To further compare experiments and simulations, it is useful to go beyond averages and map out the velocity distribution in the cytoplasm. The distribution of velocities of organelles in the cytoplasm is available indirectly from light scattering experiments measuring Doppler shift in algal cells *in vivo*: these suggest that the distribution is sharply peaked with only about 10 per cent variation close to the estimated myosin velocity [8]. Figure 3 shows the distributions of the velocities in our simulations, computed on a regular array of points inside the endoplasm, which may be ascribed to tracer organelles within the cell. It can be seen that, in the case where slip is allowed, our simulated distribution is sharply peaked, and the standard deviation is about 10 per cent or less of the mean value, in good semiquantitative agreement with experiments. In the no-slip case, the distribution of the endoplasmic velocities is still sharply peaked, but at a much lower velocity than that of the motors, at odds with experiments.

### 3.2. The flow inside the endoplasm is quasi-one dimensional

Our simulations consider a fully three-dimensional flow, with a typical simulation volume being equal to 25 × 20 × 100 in LU (see §2.3 for a mapping to real units). As we directly model the two-dimensional carpet of moving motors, together with the interaction with the wall, it is reasonable to expect that close to the motors there will be significant inhomogeneities of the fluid velocity. This may be further increased by any hydrodynamic clustering of the motors. It is then useful to ask how deep inside the endoplasm the correlations and fluctuations in the fluid velocity imparted on the intracellular solvent by the motor dynamics persist, and what the decay of such correlations is. Figure 4 shows a cut of the fluid velocity profile on different planes, stacked along the velocity gradient direction.

As can be seen in figure 4, close to the particles the flow is indeed inhomogeneous, but surprisingly weakly so. The square root of the velocity variance is only about 1–10% of the typical fluid velocity even close to the motors; it then decays exponentially with distance and becomes negligible beyond 10 LU away from the walls. In other words, our simulations show that the flow is quasi-one dimensional, corroborating the analysis of Nothnagel & Webb [11], who assumed this geometry at the outset, and also justifying the use of uniform velocity boundary conditions for the vacuole in the study of van de Meent *et al.* [5].

### 3.3. Experiments constrain the model values of viscosities, motor density and slip layer thickness

While figures 2 and 3 show that the key features of cytoplasmic streaming in plant cells may be reproduced by a reasonable set of values for cytoplasm viscosities, motor geometry, density and slip layer thickness, it should be noted that some of these control parameters are not known with high accuracy from experiments. It is therefore important to assess how our results and conclusions are affected if these are changed, within a biologically relevant range. In this section, we therefore analyse the impact of parameter variation on the fluid velocity within the endoplasm.

We first consider the effect of the thickness of the slip layer, *ɛ*. Figure 2*b* refers to *ɛ* = 220 nm, and figure 2*a* to *ɛ* = 0, with only the former data accounting for a realistic streaming within the cytoplasm. Figure 5 shows the effect of varying the slip layer, from 50 nm, which is smaller than the size of the motor, to 500 nm, the size of the cargoes. Within this slip layer, we assume that the cytosol is basically an aqueous fluid—motivated by the view that the larger viscosity inside is essentially due to macromolecular crowding with crowding agents that do not reach up to the wall, thus resulting in a lower viscosity close to the wall. The mean velocity shows a saturation behaviour as a function of *ɛ*, with essentially any value larger than 200 nm enough to account for experimental values. The lower values for slip layer thickness, while geometrically consistent, lead to a velocity of the fluid that is about 80 per cent of the motor velocity—this is slightly too low given current measurements, although the corresponding values of *ɛ* cannot be definitely discarded.

Figure 6 shows the effect of changing the density of motors, measured by their surface coverage *ϕ* = *A*_{occ}/*A*_{tot}, where *A*_{occ} is the surface area occupied by cargo beads and *A*_{tot} is the maximum area available to them without blocking cargoes on adjacent lanes (i.e. cargoes placed on a rectangular lattice just touching their neighbouring cargoes on the same lane and on adjacent lanes). It is interesting to note that, once the slip layer thickness is large enough, even a surface density equal to 0.1, i.e. 10 times smaller than full coverage, is enough to account for the experimental velocities. Figure 6 also shows that the distribution of particles on the surface plays a role. For instance, one may achieve the same surface coverage *ϕ* = 0.1 by using fewer actin lanes with more myosin motors on each. The largest ratio between fluid and particle velocity is obtained when actin lanes are closely spaced. This effect is due to the creation of channels of slower fluid between actin lanes if the lanes are spaced further apart.

Finally, figures 7 and 8 show the effects of cytoplasm viscosity and thickness, respectively. The viscosity of the cytoplasm has been measured by several authors and values in the range 6–250 cP have been reported, whereas the viscosity of the vacuole may be taken to be close to that of water, 1 cP. The cytoplasm viscosity may actually depend on probe size and this makes it difficult to assess. For our case, as organelles were used in the original experiments to track the flow, it may be more appropriate to use viscosities observed with probes of about the same size, hence towards the top end of the range quoted earlier. Figure 7 shows that the cytoplasm viscosity does have a sizeable effect on the average cytoplasmic streaming velocity: the threshold beyond which realistic streaming may be achieved via cargo dragging is about 20 cP, for which the average fluid velocity is approximately 90 per cent of the average motor velocity. Figure 8, on the other hand, shows that the thickness of the cytoplasm plays a comparatively much less important role.

## 4. Discussion and conclusions

The work we have presented provides a microscopic simulation of cytoplasmic streaming in a geometry relevant to that of plant cells, where there is an endoplasmic layer surrounding a much larger vacuole, and the flow inside the endoplasm is driven by molecular motors that run on actin lanes dragging vesicles. Therefore, our calculations can probe the viability of current continuum fluid dynamic theories of cytoplasmic streaming, which often lump the microscopic detail of the motor motion into a shear velocity boundary condition for a purely continuum problem inside the vacuole [5].

Our main conclusion is that, if motors on the cytoplasm are attached to vesicles or organelles, the drag exerted on the fluid through their motion is enough to lead to an essentially continuum flow within the endoplasm and the vacuole, with average fluid velocity fully consistent with those observed in cytoplasmic streaming. This conclusion holds for a range of realistic parameters, and suggests that there is no need to invoke a solid-like coupling between the motors moving along the subcortical actin tracks and the cytoskeletal network, or endoplasmic reticulum. Such an elastic coupling was deemed necessary in order to account for the experimental velocities in the previous work [11], and it is therefore useful to point out some key differences between our approaches. Most important is the fact that our work presents direct hydrodynamic simulations of the motion of solid spheres close to a boundary in a two-layered fluid mimicking the endoplasm–vacuole of an algal cell. As such, we fully take into account sphere–sphere and sphere–wall interactions, and crucially also consider slip at the wall. The hydrodynamic treatment in the study of Nothnagel & Webb [11], on the other hand, neglected these interactions, and this is probably at the origin of the quantitative discrepancy between our conclusions (although it should be noted that even the results of Nothnagel & Webb [11] did not definitely rule out the purely viscous coupling advocated in our work). Indeed, many-body hydrodynamics leads to a collective effect, which increases the drag exerted on the fluid by a moving carpet of motors, whereas the wall, which could counteract that effect in principle, does not counteract it in our simulations owing to the inclusion of a slip boundary condition. Our approach of introducing slip boundary conditions can probably be generalized to a range of biological solvents whose high effective viscosity is due to macromolecular crowding, so that a layer of lower effective viscosity forms close to walls, where the crowding agents (which are macromolecules and biopolymers) are depleted [20].

It is also important to recognize some remaining limitations of our calculation. First, the interior of the plant cell we wish to model is actually cylindrical, while we have modelled a slab of fluid. This is perhaps not too important as the radius of curvature of the cell is much larger than the thickness of the endoplasm, which is the region of interest, as this is where the coupling between motor motion and vacuolar fluid dynamics occurs. Second, we have neglected all details of the intracellular fluid and just modelled the cytoplasm/endoplasm as a Newtonian viscous fluid. In reality, the cytoplasm is non-Newtonian, shear thinning, and it would be interesting to see how a more careful model of the intracellular solvent may affect our results. To this end, we should however make a decision as to whether to model the endoplasm as simply a crowded medium or to include the endoplasmic reticulum as a polymer network. Given the uncertainty with which several of the parameters we use are known, this does not yet seem appropriate and might have to await more detailed experimental study of the rheological propertites of the endoplasm. In this context, it is also interesting to note that in the work of Niwayama *et al.* [4], who studied cytoplasmic streaming in *C. elegans*, simulations at constant viscosity proved sufficient to quantitatively reproduce the streaming motion of the (actually non-Newtonian) cytoplasm. A final simplification of our treatment is that the membrane separating endoplasm and vacuole was computationally treated as an infinitesimally thin interface where the fluid velocity was continuous, whereas, in reality, this membrane is better described as a viscoelastic membrane.

## Acknowledgements

This work was supported by a DAAD post-doctoral fellowship to K.W. and by EPSRC grant EP/E030173. M.E.C. holds a Royal Society Professorship.

## Appendix A

In this appendix, we give additional information on our numerical framework and on the derivation of the formulae determining the partial slip parameter in our simulations (see equations (2.5)–(2.7)).

In order to model slip at a solid wall in LB simulations, a partial slip parameter *p* can be introduced that specifies the reflectivity of the wall. This parameter determines the fraction of a velocity population that is reflected specularly (or bounced forward) at the wall in an LB streaming step (figure 9). Note that *p* = 0 corresponds to full no-slip boundary conditions for the velocity field at the wall, whereas *p* = 1 corresponds to full slip—we refer to any value of *p* in between as *partial slip*. The wall reflectivity can be interpreted physically in terms of the Knudsen number, which enters the simulation via the relaxation time *τ* [17]. Another possibility is to interpret the wall slip in a more general way where *p* is simply an open parameter that can be determined if the slip length *u*(0)/*u*′(0) is known from experiment or theory.

One particular setting that can be described in this way is that of a two-viscosity system where a thin low-viscosity layer is situated close to the wall and a bulk fluid of higher viscosity beyond. If the behaviour of the thin layer itself is not of interest, all the layer does is introduce effective slip for the bulk system which can be simulated using the slip parameter *p*.

Assuming the velocity profile of the fluid to be of the form
the slip length *u*(0)/*u*′(0) can be determined analytically for *ɛ* ≪ *R*, where *R* is the total system size and *ɛ* is the thickness of the low-viscosity layer. Various conditions on *u* give
and therefore

Observing that *u*(0) = *b*_{0} and *u*′(0) = *b*_{1}, we arrive at
A1to first order in *ɛ*, where *η*_{I} is the (lower) viscosity of the thin layer close to the wall and *η*_{II} is the (higher) viscosity of the bulk system.

In order to relate slip length *u*(0)/*u*′(0) and slip parameter *p* and relaxation time *τ* in the LB method, we consider the D2Q9 model. All simulations have been performed in the D3Q19 model where the derivation works in exactly the same way but the algebra is more cumbersome (figure 10).

In a collision step, velocity populations of the same node collide and relax towards their local and instantaneous equilibrium population *f*_{i}^{eq} (*t*, **x**)
A2

Here, *f*_{i}(*t*, **x**) is the velocity population of direction *i* at time *t* and node **x** after collision and *f*_{i}^{*}(*t*, **x**) is the population before collision. The relaxation time *τ* determines how far the population relaxes in each step and is related to viscosity by *τ* = (1 + 6*η*_{II})/2. Here only the bulk viscosity *η*_{II} enters the formula as it is the bulk fluid that is simulated in LB, whereas the low-viscosity layer, and hence *η*_{I}, is incorporated into the boundary condition.

The equilibrium populations are given by A3

with *ω*_{1} = *ω*_{3} = *ω*_{5} = *ω*_{7} = 1/9 = :*w*_{1} and *ω*_{2} = *ω*_{4} = *ω*_{6} = *ω*_{8} = 1/36 = :*w*_{2}. Density *ρ* is assumed to be constant, and fluid velocity is given by
A4

We next need to make a couple of assumptions on the velocity profile. First, we assume that **u**(*t*, **x**) = *u*(*t*,*z*)**e**_{x} is parallel to the wall and that it only depends on the distance to the wall *z*. So, for the velocity in the *x*-direction, we have
A5

In the streaming step, velocity populations from time step *t* are passed on to neighbouring nodes and become pre-collision populations of time step *t* + 1. Below, streaming is given for the layer of fluid nodes just above the wall (i.e. those experiencing reflection from the wall as follows). As velocities are assumed to be homogeneous in the *x*- and *y*-directions, those indices are omitted altogether and *f*_{i,k} denotes velocity population *i* in layer *k* above the first fluid layer. For the first layer, *k* = 0 will be omitted,
A6

The next assumption to be made is that a steady state has been reached and the velocity profile and all populations no longer change with time, meaning that all reference to *t* can be dropped.

Plugging streaming (equation (A 6)) and collision (equation (A 2)) into equation (A 5) yields

To close this equation, we need to express the pre-collision populations *f*_{i}^{*} in terms of equilibrium populations *f*_{i}^{eq} as for those we have equation (A 3) relating them to *ρ* and *u*.

For populations *f*_{1}^{*} and *f*_{5}^{*}, this is easy as no populations from higher layers (larger *z*) are involved and the profile is assumed to be homogeneous in *x*, meaning that after many iterations *f*_{1}^{*} and *f*_{5}^{*} simply relax to *f*_{1}^{eq} and *f*_{5}^{eq}. For *f*_{6}^{*} and *f*_{8}^{*}, this is more complicated, as populations from higher and higher layers will be involved
A7

For large enough *n*, the term involving can be disregarded and using the earlier mentioned expression and equation (A 3) in the equation for *ρ**u* leaves us with
where *u*_{0} is the velocity in the layer of fluid nodes closest to the wall and *u*_{k} is the velocities in higher layers.

Using Taylor expansions for velocities at higher layers and combining subsequent terms of *u*_{k} and *u*_{k+1} in the sum above finally results in
A8

In principle, can be evaluated for arbitrary *j* but then more and more derivatives *u*_{0}^{(j)} have to be known in advance. We therefore truncated after the linear term and checked that the resulting formula still holds within our numerical precision for Poiseuille flow where the velocity profile in the channel is quadratic.

In the linear approximation *u*_{0}^{(j)} = 0 for *j* > 1, we obtain

At this point, it is important to note that, in line with standard LB algorithms working within the Euler scheme that we use, the wall is implemented between two nodes, so that it is mid-grid. Therefore, *u*(0) from equation (2.5) does not correspond to the velocity at the first layer of fluid nodes *u*_{0} but lies half a lattice spacing below *u*(0) = *u*_{0} − 1/2*u*′_{0}. As *u*′_{0} is assumed constant, it is of course the same as *u*′(0) and we get
A9which can be used together with equation (A 1) to obtain the formula relating slip parameter *p* and physical quantities *ɛ*, *η*_{I} and *η*_{II} in equation (2.7) of §2.

## Footnotes

↵1 Cytoplasmic streaming may also be observed in non-plant cells—for instance, in the nematode

*Caenorhabditis elegans*[4].↵2 This is typical in LB work [15,18]. Because the LB algorithm uses fluid inertia to update the dynamics, the most efficient discretization uses the highest density compatible with still remaining in the low Reynolds number regime. With our chosen density, the Reynolds number associated with the vesicle size never exceeds 0.1.

- Received December 12, 2011.
- Accepted January 26, 2012.

- This journal is © 2012 The Royal Society