Enhanced Transport of Two Spheres in Viscous Fluid

TL;DR

This paper presents a numerical study showing that two spheres moving in viscous fluid can travel farther together than alone, implying hydrodynamic interactions can enhance transport efficiency in biological systems.

Contribution

It introduces a numerical solution for the synchronized motion of two spheres in viscous fluid, revealing how their interaction improves transport efficiency.

Findings

Presence of a second sphere increases travel distance for a given work

Hydrodynamic interactions can enhance molecular motor cargo transport

Results suggest potential biological implications for cargo movement

Abstract

We obtain a numerical solution for the synchronous motion of two spheres moving in viscous fluid. We find that for a given amount of work performed, the final distance travelled by each sphere is increased by the presence of the other sphere. The result suggests that the transport efficiency of molecular motor cargo in vivo may be improved due to an effective hydrodynamic interaction with neighboring cargos moving along the same direction.

Equations

$$F_{\mathrm{fluid}}$$

$$\left({\rho }_s+\frac{{\rho }_f}{2(1+3{ϵ}^3)}\right)\dot{v}=-\frac{9\eta }{2R^2}\frac{v}{1+\frac{3}{2}ϵ-{ϵ}^3}$$

$$-\frac{9\eta }{2R^2}{\int }_{-\infty }^{t}d{t}^{\prime }\dot{v}(t^{\prime })h_{ϵ}(t-t^{\prime })+\frac{3}{4\pi R^3}F(t),$$

$$({\rho }_s+\frac{{\rho }_f}{2})\dot{v}=-\frac{9\eta }{2R^2}v$$

$$-\frac{9}{2R}\sqrt{\frac{\eta {\rho }_f}{\pi }}{\int }_{-\infty }^{t}d{t}^{\prime }\frac{\dot{v}(t^{\prime })}{\sqrt{t-t^{\prime }}}+\frac{3}{4\pi R^3}F(t),$$

$$\theta$$

$$f(\theta )$$

$$\frac{2{\rho }_s(1+3{ϵ}^3)+{\rho }_f}{(1+3{ϵ}^3)(2{\rho }_s+{\rho }_f)}\frac{du}{d\theta }=-\frac{u}{1+\frac{3}{2}ϵ-{ϵ}^3}$$

$$x(\theta )$$

$$w(\theta )$$

$$f_{\mathrm{drive}}\equiv \frac{w(\theta )}{\Delta x(\theta )}$$

$$z(\theta )\equiv \frac{f_{\mathrm{drive}}(\theta )}{\bar{u}(\theta )}=\frac{w(\theta )\Delta \theta }{\Delta x(\theta {)}^2}.$$

$$f_{\rm drive}(\theta)\equiv\frac{w(\theta)}{x(\theta)}=\left\{\begin{array}[]{l}1,\ (\theta<\theta_{\rm pulse})\\ x(\theta_{\rm pulse})/x(\theta),\ (\theta\geq\theta_{\rm pulse}).\end{array}\right.$$

$$T_{\mathrm{pulse}}=10\mu \mathrm{s}=20{\tau }_B=20\times \frac{3\times 10^3\mathrm{kg}/{\mathrm{m}}^3\times R^2}{9\times 2\times 10^{-3}kg/(\mathrm{m}\cdot \mathrm{s})},$$

$$R\simeq 2\mu \mathrm{m}\sim O(1\mu \mathrm{m}),$$

$$q\equiv \frac{Q(t)⟨\delta X(t{)}^2⟩}{⟨X(t){⟩}^2}$$

$$h_{ϵ}(t)=\frac{1}{\pi }{\int }_0^{\infty }{\hat{h}}_{ϵ}(s)e^{-st/{\tau }_{\nu }}ds$$

$${\hat{h}}_{ϵ}(s)=\mathrm{Im}\left[\frac{A(s)}{B_{ϵ}(s)}\right]$$

$$A(s)$$

$$B_{ϵ}(s)$$

$$h_{ϵ}(t)$$

$$F_{\mathrm{fluid}}$$

$$g_{ϵ}(t)=\frac{1}{\pi }{\int }_0^{\infty }{\hat{g}}_{ϵ}(s)e^{-st/{\tau }_{\nu }}ds$$

$${\hat{g}}_{ϵ}(s)=\mathrm{Im}\left[\frac{A(s)}{C_{ϵ}(s)}\right]$$

$$A(s)$$

$$C_{ϵ}(s)$$

$$g_{ϵ}(t)$$

$$a(\theta )=-F(ϵ)G(ϵ)u(\theta )-F(ϵ){\int }_0^{\theta }d{\theta }^{\prime }a({\theta }^{\prime }){\tilde{h}}_{ϵ}(\theta -{\theta }^{\prime })+F(ϵ)f(\theta ).$$

$$a(\theta )$$

$${\tilde{h}}_{ϵ}(\theta )$$

$$F(ϵ)$$

$$G(ϵ)$$

$${\int }_{a+k\Delta x}^{a+(k+1)\Delta x}f(x)dx\simeq \Delta x\left[\frac{1}{2}f(a+k\Delta x)+\frac{1}{2}f(a+(k+1)\Delta x\right]$$

$$f(x)\sim \frac{A}{(b-x{)}^{\alpha }}(0<\alpha <1)$$

$${\int }_{b-\Delta x}^{b}f(x)dx\simeq {\int }_{b-\Delta x}^b\frac{A}{(b-x{)}^{\alpha }}=\frac{A(b-\Delta x{)}^{1-\alpha }}{1-\alpha }.$$

Full text

Subcellular transport processes have been the subject of continuing interest largely because of the high efficiency of molecular-scale motors responsible for such transport, as compared to the efficiency of their macroscopic counterparts Pietzonka et al. (2016a); Jr et al. (2000); Wang and Oster (2002a); Schmiedl and Seifert (2008a); Lau et al. (2007); Ariga et al. (2018); Kolomeisky (2013); Kolomeisky and Fisher (2007); Neuman et al. (2005); Hwang and Hyeon (2017); Bustamante et al. (2001); Wang and Oster (2002b); Schmiedl and Seifert (2008b); Brown and Sivak (2017); Hwang and Hyeon (2018). An important factor affecting such a transport process is the hydrodynamic interactions Zhu et al. (1992); Kao et al. (1993); Henderson et al. (2002); Liverpool and MacKintosh (2005); Atakhorrami et al. (2005); Kheifets et al. (2014); Alder and Wainwright (1970); Kim and Matta (1973); Fedele and Kim (1980); Parmar et al. (2011); Lesnicki and Vuilleumier (2016); Jung and Schmid (2017); Arminski and Weinbaum (1979); Ardekani and Rangel (2006); Marchetti et al. (2013); Salbreux et al. (2009); Padding and Louis (2004); Wu et al. (2018); Cattuto et al. (2006); Detcheverry and Bocquet (2012) of the molecular motor and cargo with the surrounding viscous fluid Goldtzvik et al. (2016). First, it has been shown that the Stokes drag experienced by a spherical object embedded in a fluid is reduced in the presence of another sphere in the vicinity Stimson and Jeffery (1926). The reduction of the drag can be understood intuitively as coming from the indirect transfer of the momentum between the embedded objects, mediated by the fluid. Secondly, there is a separate effect of hydrodynamic memory due to the unsteady flow Boussineq (1885); Basset (1888); Oseen (1927); Kheifets et al. (2014); Alder and Wainwright (1970); Kim and Matta (1973); Fedele and Kim (1980); Parmar et al. (2011); Lesnicki and Vuilleumier (2016); Jung and Schmid (2017); Arminski and Weinbaum (1979); Ardekani and Rangel (2006), where the momentum transferred from an embedded object is transiently stored in the fluid, and then transferred back to the object at a later point Arminski and Weinbaum (1979). Such an effect may facilitate the transport even for a single spherical object Seyler and Pressé .

These considerations motivate us to ask whether fluid flow generated by the motion of a moving object in viscous fluid indeed facilitates the transport of a neighboring cargo. To address this question here, we examine the displacement of two neighboring spheres under the same constant external force of finite duration by numerically integrating the equation of motion Ardekani and Rangel (2006) (Eq.(2)).

We will show that for a given distance traveled, the required input work for two neighboring spheres is less than that for spheres separated by large distances. We also find that two neighboring spheres are transported faster than separated spheres, supporting the idea that indirect, fluid-induced, inter-vesicle interaction may be utilized in active subcellular transport for improved efficiency.

Equation for two spheres– The equation of motion for two spheres moving in viscous fluid has previously been derived by generalizing the equation for one sphere in a fluid Ardekani and Rangel (2006). Briefly, fluid flow is computed in the presence of a sphere of radius $R$ with no-slip boundary condition by using the unsteady Stokes equation where the convective term in the Navier-Stokes equation is neglected but the time-derivative is kept Landau and Lifshitz (1987). Another sphere is introduced at a center-to-center distance $d$, then the flow is perturbed in order to satisfy the no-slip boundary at the second sphere, and perturbed again in order to satisfy the boundary condition at the first sphere. This process is iteratively repeated. The series is truncated so that the resulting error is of order $\epsilon^{4}$ where $\epsilon\equiv R/d$ Ardekani and Rangel (2006). The resulting expression for the force has a remarkable symmetry with respect to the exchange of the spheres’ positions. In particular, when their radii are the same and they move with the same velocity ${\bf v}(t)$, the force exerted on each of the spheres is the same. In the special case of one-dimensional motion, along a line connecting their centers, with the same velocity component $v(t)$ (Fig. 1), the equation reduces to Ardekani and Rangel (2006)

[TABLE]

where $\rho_{f}$ and $\eta$ are the density and the dynamic viscosity of the fluid, respectively. The functional form of the memory kernel $h_{\epsilon}(t)$ is formally given as a Laplace transform Ardekani and Rangel (2006); see SI for details.

A similar expression can be obtained for the spheres whose motion is in the direction perpendicular to the lines joining their centers (see SI). Because the force exerted on each sphere by the fluid is the same, the motion of both spheres can be kept synchronous by applying the same external force $F(t)$ under the same initial velocities. The equation of motion for each sphere can be written as

[TABLE]

where $\rho_{s}$ is the density of the sphere. By comparing with the exact result Stimson and Jeffery (1926) for the special case of the motion with constant velocity, it was argued that Eq. (2) is a reasonable approximation for $\epsilon\lesssim 0.25$  Ardekani and Rangel (2006). The functional form of $h_{\epsilon}(t)$ is such that $h_{0}(t)=R\sqrt{\frac{\rho_{f}}{\pi\eta t}}$ (see SI), so that for $\epsilon=0$, we obtain the familiar Basset-Boussinesq-Oseen equation Boussineq (1885); Basset (1888); Oseen (1927); Kheifets et al. (2014) for an accelerating sphere in a fluid

[TABLE]

where the second term on the right hand side captures the effect of the hydrodynamic memory Boussineq (1885); Basset (1888); Oseen (1927); Kheifets et al. (2014).

For convenience, we may rewrite the equation for synchronous movement of two spheres in terms of dimensionless quantities, defined as

[TABLE]

where $F_{\rm max}$ is the maximum value of $F(t)$ and $\tau_{B}\equiv(2\rho_{s}+\rho_{f})R^{2}/9\eta$ is the Brownian relaxation time Padding and Louis (2006). Then, Eq. (2) is rewritten as

[TABLE]

The velocity $u(\theta)$ can be obtained by numerical integration of Eq. (5), by discretizing the time $\theta$. The integral can then be performed by a simple trapezoidal rule, which gives a reasonably accurate result when the discretization is performed with step $\Delta\theta=0.001$, as can be checked for the case of $\epsilon=0$ where exact solutions are available for certain special cases (see SI). Once $u(\theta)$ is obtained, it is then straightforward to obtain the non-dimensionalized quantities $x$ and $w$ corresponding to position and work, respectively, by additional integration:

[TABLE]

where $X(t)$ is the position of the sphere. In order to quantify a notion of efficiency for transport, we define the dimensionless effective transport drive force $f_{\rm drive}$ Arminski and Weinbaum (1979),

[TABLE]

where $\Delta x(\theta)\equiv x(\theta)-x(0)$ is the displacement. This quantity is the dimensionless version of the specific energy consumption, often used for the measure of transport efficiency in the transportation industry Gabrielli and von Kármán (1950); Trancossi (2016); Chiara et al. (2017). The smaller value of $f_{\rm drive}(\theta)$ implies less amount of external work required for a given displacement. We also define the dimensionless effective friction Arminski and Weinbaum (1979)

[TABLE]

Improved transport of two neighboring spheres– We now consider a simple protocol where a constant external force is applied over a finite duration $T_{\rm pulse}$, starting from $t=0$. Note that by definition, the maximum value of the normalized force $f(\theta)$ is unity. Therefore, $f(\theta)=1$ for $0\leq\theta\leq\theta_{\rm pulse}$ and zero otherwise, where $\theta_{\rm pulse}\equiv T_{\rm pulse}/\tau_{B}$. We also take $x(0)=u(0)=0$. Therefore, the input work is $w(\theta)=f\times x(\theta)=x(\theta)$ for $\theta<\theta_{\rm pulse}$ and $w=x(\theta_{\rm pulse})$ otherwise. As such,

[TABLE]

The results of the numerical computation for $\theta_{\rm pulse}=20$ and $\rho_{s}=\rho_{f}$ are shown in Fig. 2, where the values of $a\equiv\frac{du}{d\theta}$, $u$ and $x$ are compared for $d=\infty$, $d=8R$, and $d=4R$. The case of $d=\infty$ corresponds to the motion of a single sphere. We find that for a given strength of the external force, the magnitudes of both acceleration and deceleration for two neighboring spheres are larger than those for the spheres separated with a larger distance, and the overall effect is such that $v(t)$ for neighboring spheres is larger for all values of $t$.

More importantly, $x(\infty)$ increases by a large amount when the inter-sphere distance $d$ decreases, whereas $x(\theta_{\rm pulse})$ is only weakly dependent on distance. Therefore, from the second line of Eq. (9), we find that $f_{\rm drive}(\infty)$ for neighboring spheres is less than that of a separate sphere. That is, the neighboring spheres travel farther compared to the separate spheres, for a given amount of input work. The graphs of $f_{\rm drive}(\theta)$ are shown for several values of $d$ in Fig. 3, where we see that in fact $f_{\rm drive}(\theta)$ is an increasing function $d$ for all values of $\theta$. The values of $f_{\rm drive}(\infty)$ are also plotted in Fig. 4 for several values of $R/d$ and $\theta_{\rm pulse}$, where approximate values of $x(\infty)$ are obtained from the values of $x(\theta)$ in those approximately flat regions of the $x$ plots (Fig. 2) with $u(\theta)\leq 0.01$. The trend of the reduced effective drive force for neighboring spheres is evident. In particular, we see that for $\theta_{\rm pulse}=20$, the value of $f_{\rm drive}$ reduces from 0.682 at $R=\infty$ to 0.596 at $d=4R$, resulting in about $\sim 13$ % reduction in the required work for a given displacement. The reduction of $f_{\rm drive}$ is smaller for larger values of $\theta_{\rm pulse}$. This is because the motion of a sphere is diffusive at the time scale of $t\gg\tau_{B}$ Kheifets et al. (2014), with only tiny effects coming from hydrodynamic memory. This can be seen from the graphs of $x(\theta)$ for $\theta_{\rm pulse}=1000$, shown in the inset of Fig. 2 (c), where the slope is approximately proportional to the instantaneous applied force, showing the typical behavior of an overdamped particle. In particular, since the motion of the sphere almost stops after the force pulse, we have $x(\theta_{\rm pulse})\simeq x(\infty)$, leading to $f_{\rm drive}(\infty)=x(\theta_{\rm pulse})/x(\infty)\simeq 1$ regardless of the inter-sphere distance, implying negligible reduction of the effective transport drive force.

We also plot the graph of $z(\theta)$ in Fig. 5 for $d=\infty$ (single sphere), $d=8R$, and $d=4R$, showing that for given displacements, two neighboring spheres not only requires less work but also results in faster transport as compared to separate spheres. The value of $z(\theta)$ here diverges as $\theta\to\infty$ because $\lim_{\theta\to\infty}\bar{u}(\theta)=0$, but it will be maintained at finite values when a periodic force is applied so that a nonequilibrium steady state is reached Arminski and Weinbaum (1979).

Order of magnitude estimates with biological parameters– To see whether the reduced effective transport drive force and friction due to the indirect interaction between two spheres is relevant for subcellular transport processes, we perform an order of magnitude estimate using biological parameters. More specifically, we consider an example of cargo transport by a kinesin motor, where the constant force pulse can be considered as an extremely simplified model of the force exerted by the kinesin motor and its cargo. The force duration may be taken as $T_{\rm pulse}=10\ {\rm\mu s}$, the time scale during which the stepping motion occurs Carter and Cross (2005); Zhang and Thirumalai (2012); Goldtzvik et al. (2016). As discussed previously, the hydrodynamic memory plays a role only if $T_{\rm pulse}$ is not too much larger than $\tau_{B}\equiv(2\rho_{s}+\rho_{f})R^{2}/9\eta$. For fixed values of $T_{\rm pulse},\rho_{s},\rho_{f}$, and $\eta$, this tells us that the size of the cargo should be sufficiently large in order for hydrodynamic memory to play a significant role. For example, we previously found that $f_{\rm drive}(\infty)$ for two spheres with $d=4R$ was 13 % less than for those separated by an infinite distance. Using the values $\eta=2\times 10^{-3}\ {\rm kg/(m\cdot s)}$ Ando and Skolnick (2010); Luby-Phelps (2000); Ridgway et al. (2008); Verkma (2002) (BNID150903)111The ID number of BioNumbers Database Milo et al. (2010). and $\rho_{s}=\rho_{f}=10^{3}\ {\rm kg}/{\rm m}^{3}$ Heyden and Ortiz (2017); Moran et al. (2010) (BNID113851) for the cytoplasm, we get

[TABLE]

leading to

[TABLE]

about the size of a large size vesicle such as an organelle Lodish et al. (2000); Hirokawa (1998); Hancock (2014). Considering the fact that organelles of size of order $1\ \mu\rm m$ are often transported by molecular motors Lodish et al. (2000); Hirokawa (1998); Hancock (2014), the reduction of effective drive force and friction driven by hydrodynamic effects warrants further investigation within living environments.

Discussion– In this work, we presented a numerical solution of two spheres moving in synchrony in a viscous fluid with the same force applied to each sphere. We found that for a given displacement for each sphere, the required work is less for two neighboring spheres than for spheres separated by a large distance, and the former is transported faster than the latter. In reality, the asynchrony in cargo transport may somewhat reduce the effect proposed here. Study of such a generalized case is straightforward, albeit technically more involved.

Our results support the idea that the efficiency of subcellular transport may be improved by hydrodynamic interaction between the neighboring cargos. Taking thermal fluctuations explicitly into account, the transport efficiency is quantified by  Hwang and Hyeon (2018); Dechant and Sasa (2018)

[TABLE]

where $Q(t)$ is the energy consumption up to time $t$, $\langle X(t)\rangle$ is the average displacement, and $\langle\delta X(t)^{2}\rangle$ the transport precision. It has been shown that $q\geq 2k_{B}T$, and this fundamental bound is called the thermodynamic uncertainty principle Barato and Seifert (2015); Gingrich et al. (2016); Pietzonka et al. (2016b); Pigolotti et al. (2017); Hyeon and Hwang (2017); Proesmans and Van den Broeck (2017). Within this bound, the molecular motor that performs transport with minimal energy expenditure and with highest precision is the most efficient one by definition. Since we considered the solution to a deterministic equation, the displacement we obtained is expected to be the thermally averaged displacement. Since the input work is proportional to $Q(t)$ for a given value of thermal efficiency, and since our results tell us that $\langle X\rangle$ for two neighboring spheres is larger than that for separate spheres, Eq.(12) tells us that two neighboring spheres have higher transport efficiency due to hydrodynamic interactions, if $\langle\delta X(t)^{2}\rangle$ are the same for both cases. Full analysis of the transport efficiency, taking into account transport precision in the presence of the thermal fluctuations, would require more sophisticated formalism such as the fluctuating hydrodynamics Padding and Louis (2004); Wu et al. (2018); Cattuto et al. (2006); Detcheverry and Bocquet (2012).

I Acknowledgement

JL was supported by the National Research Foundation of Korea, funded by the Ministry of Education (NRF-2017R1D1A1B03031344). SS and SP were supported by ARO grant W911NF-17-1-0162 on “Multi-Dimensional and Dissipative Dynamical Systems: Maximum Entropy as a Principle for Modeling Dynamics and Emergent Phenomena in Complex Systems”.

Appendix A Memory kernel for the two sphere equation (Eq. (2))

The memory kernel in Eq. (5) is given by Ardekani and Rangel (2006)

[TABLE]

where

[TABLE]

with

[TABLE]

and $\tau_{\nu}\equiv\rho_{f}R^{2}/\eta$. We note that for $\epsilon\to 0$, contribution from the term with the factor $\exp\left(-i\sqrt{s}\left(1-\epsilon^{-1}\right)\right)$ vanishes in the integral Eq. (13), and we have

[TABLE]

recovering the memory kernel for one sphere Boussineq (1885); Basset (1888); Oseen (1927).

Appendix B The equation for spheres moving perpendicular to their line of centers

In the main text, we focused on the case where the spheres move along the line connecting their centers. We may also consider two identical spheres moving perpendicular to the center-to-center axis, as shown in Fig. 6. The resulting force exerted by the fluid on each sphere is Ardekani and Rangel (2006)

[TABLE]

where $g_{\epsilon}(t)$ is formally given as a Laplace transform Ardekani and Rangel (2006),

[TABLE]

where

[TABLE]

with

[TABLE]

and $\tau_{\nu}\equiv\rho_{f}R^{2}/\eta$. We note that for $\epsilon\to 0$, contribution from the term with the factor $\exp\left(-i\sqrt{s}\left(1-\epsilon^{-1}\right)\right)$ vanishes in the integral Eq. (13), and we have

[TABLE]

recovering the memory kernel for one sphere Boussineq (1885); Basset (1888); Oseen (1927). The qualitative behavior for two neighboring spheres moving in the direction perpendicular to the central line is similar to those moving along the central line, in that the transport is enhanced, as shown in Figs. 7 and 8, where bth positions and the transport drive force are compared. We again find about 13 % reduction in transport drive force for $d=4R$.

Appendix C Solution of integro-differential equation

We numerically solved the integro-differential Eq. (5), rewritten as

[TABLE]

where we assumed that $v=0$ for $\theta\leq 0$, with

[TABLE]

The integral in Eq. (22) is performed by using the trapezoidal rule, where we compute an integral of the form $\int_{a}^{b}f(x)dx$ by first discretizing the interval $[a,b]$ into $N$ subintervals, and approximating the integral in each subinterval as an area of the trapezoid,

[TABLE]

where $\Delta x\equiv(b-a)/N$. Care must be taken when the function $f(x)$ is divergent at a boundary, say at $x=b$. Here, we cannot use Eq. (24) at the subinterval $[b-\Delta x,b]$, so the corresponding integral must be treated separately. We utilize the asymptotic form for the function $f(x)$ for $x\to b$,

[TABLE]

to make the approximation

[TABLE]

By using Eqs. (24) and (26), we obtain

[TABLE]

For the integral in Eq. (22), we have

[TABLE]

The last integral must be treated separately using Eq. (26) since $\tilde{h}_{\epsilon}(\theta)$ diverges at $\theta=0$. Therefore we have to find the asymptotic form $\tilde{h}_{\epsilon}(\theta)$ for $\theta\to 0$. In fact, from Eq. (13), we see that $h_{\epsilon}(t)$ is dominated by the value of $\hat{h}_{\epsilon}(s)$ at $s\to\infty$ as $t\to 0$. The oscillatory contribution containing the factor $e^{-i\sqrt{s}(1-\epsilon^{-1})}$ vanishes as $s\to\infty$, and we have

[TABLE]

and consequently

[TABLE]

Note that the right-hand side of Eq. (22) contains $u(\theta)$ that is undetermined at the time when computing $a(\theta)$. It is to be computed using the trapezoidal rule

[TABLE]

We simply substitute Eq. (31) into Eq. (22), along with and Eqs. (28) and (30) to get

[TABLE]

Moving the term proportional to $a(\theta)$ to the right-hand side to the left-hand side and solving for $a(\theta)$, we get

[TABLE]

where

[TABLE]

Eq. (33) allows us to compute $a(\theta)$ in terms of $u(\theta-\Delta\theta)$ and $a(0),a(\Delta\theta),\cdots a(\theta-\Delta\theta)$, by storing $a(\theta)$’s as arrays during the computation. Once $a(\theta)$ is obtained, $u(\theta)$ can be obtained by Eq. (31).

We also need to compute the integral in Eq. (13) in order to obtain $\tilde{h}_{\epsilon}(\theta)$. The upper limit of the integral is infinity, so we truncate the integral when the integrand is sufficiently small. In other words, we truncate the region with $e^{-s\tau_{B}\theta/\tau_{\nu}}<\delta$. The remaining integral is obtained numerically by the trapezoidal rule. It would be computationally inefficient to perform the integral in Eq. (13) each time we compute $a(\theta)$. Therefore, we compute $\tilde{h}_{\epsilon}(0),\tilde{h}_{\epsilon}(\Delta\theta),\cdots\tilde{h}_{\epsilon}(N\Delta\theta)$ at the start of the computation and store them as arrays, where $N\Delta\theta$ is the upper limit of $\theta$ that $a(\theta)$ will be computed.

For an isolated sphere, an analytic solution of Eq. (3) for a constant pull is known Arminski and Weinbaum (1979) which can be compared with our numerical solution for $\epsilon=0$ to assess the accuracy of our method. We found that using $\Delta\theta=0.001$, $\Delta s=0.0001$, and $\delta=0.0001$ yields reasonably accurate solution, as can be seen in Fig. 9 where the numerical and the analytic solutions for the velocities are compared for $\rho_{s}=0$222The analytic solution for $\rho_{s}=\rho_{f}$ is expressed in terms of error function with complex argument Arminski and Weinbaum (1979), For convenience, we therefore use $\rho_{s}=0$ for the purpose of plotting the analytic solution. . These parameters were also used for performing the numerical integration of the two-sphere equation.

C.1 Assessment of the accuracy of the truncation in Eq. (1)

As mentioned in the main text, Eq. (5) is obtained by a truncation where the error is considered to be of order $O(\epsilon^{4})$ where $\epsilon\equiv R/d$. It was argued in Ref. Ardekani and Rangel (2006) that Eq. (5) is reasonably accurate for $\epsilon\lesssim 0.25$, by comparing the Stokes drag in Eq. (5),

[TABLE]

with the known exact form Stimson and Jeffery (1926)333There is an error in the overall prefactor of the final expression in Ref. Stimson and Jeffery (1926) where $2/3$ is given instead of $4/3$. We could check that $4/3$ is the correct one by carefully following their derivation.

[TABLE]

where

[TABLE]

The numerical result is quite robust as we replace $f^{\rm(trunc)}_{\rm Stokes}$ by $f^{\rm(exact)}_{\rm Stokes}$ up to $\epsilon=0.25$, as shown in Fig. 10 where the results for the position $x$ are compared for $t=20\tau_{B}$ and $\rho_{s}=\rho_{f}$. This provides circumstantial evidence that Eq. (5) is quite accurate for $\epsilon\lesssim 0.25$. We however do not use $f_{\rm Stokes}^{\rm exact}$ as we avoid ad hoc mixing of exact and truncated expressions.