# Numerical and analytical approaches to an advection-diffusion problem at small Reynolds number and large Péclet number

###### Abstract

Obtaining a detailed understanding of the physical interactions between a cell and its environment often requires information about the flow of fluid surrounding the cell. Cells must be able to effectively absorb and discard material in order to survive. Strategies for nutrient acquisition and toxin disposal, which have been evolutionarily selected for their efficacy, should reflect knowledge of the physics underlying this mass transport problem. Motivated by these considerations, in this paper we discuss the results from an undergraduate research project on the advection-diffusion equation at small Reynolds number and large Péclet number. In particular, we consider the problem of mass transport for a Stokesian spherical swimmer. We approach the problem numerically and analytically through a rescaling of the concentration boundary layer. A biophysically motivated first-passage problem for the absorption of material by the swimming cell demonstrates quantitative agreement between the numerical and analytical approaches. We conclude by discussing the connections between our results and the design of smart toxin disposal systems.

## I Introduction

The survival of any microorganism is dependent on its ability to control the movement of nutrients and waste between the cell interior and the external environment. Nutrients must be easily captured and waste easily disposed with minimal metabolic exertion cole2013cost (). Factors which determine the efficiency of any material transfer mechanism for gathering and disposing of small molecules include the microscale flow profile surrounding the cell short2006flows (), the diffusion coefficient of molecules being transported, and the location of transport receptors on the cell surface whalen2012actin (); Song (). As with any other trait, the physical architecture of the transport receptor arrangement is evolutionarily selected to work efficiently within the particular environment that the organism lives. A more detailed understanding of the mass transport problem for molecules in the extracellular fluid can yield valuable information for discerning why a given architecture is preferable over other possible choices vogel1994life ().

An example which motivates the current study is understanding the mechanism of toxin disposal by sea urchin embryos whalen2012actin (); gokirmak2012localization (). Early in development, microvilli (microscopic cylindrical cell membrane protrusions) on the embryo grow to several microns in length, with transport receptors localized at the tips of the microvilli. This localization may help facilitate more efficient toxin disposal by the embryo. By releasing molecules away from the body of the cell, the chances are reduced that expelled toxins are reabsorbed. In this setting, the mass transport problem for molecules in the extracellular fluid is described an advection-diffusion equation at small Reynolds number and large Péclet number. Interestingly, the characteristic lengthscale for the concentration boundary layer may provide a physical rationale for the length of the microvilli involved in the toxin transport urchin ().

Calculating the efficiency of a material transfer mechanism requires the solution of two distinct problems. The first problem is to determine the flow profile of the extracellular fluid. In the present context, the Reynolds number of the flow is much smaller than one, which justifies the use of the Stokes equation for describing microscale flows around individual cells felderhof2015stokesian (). The second problem is to find a solution of the advection-diffusion equation using the calculated flow profile as input. The solution of the advection-diffusion equation provides the concentration profile governing the average motion of particles dissolved in the fluid. Knowledge of the concentration profile leads directly to information about the probability for a given particle to contact the outer envelope of a cell. Since the absorption of a particle requires that it first make contact with the cell surface, the concentration profile can be used to find a theoretical maximum absorption probability berg1977physics ().

Although the equations which govern the extracellular flow and diffusion of small molecules are well understood, deriving their solutions or even their approximate solutions is often quite involved ifm (); ifac (); acrivos1960solution (); acrivos1962heat (); acrivos1965asymptotic (); goddard1966asymptotic (). The difficulties faced in finding accurate analytic solutions to the advection-diffusion equation can make numerical schemes an attractive method for finding concentration profiles cfmht (); ecfd (); icfme (). In this paper we present both numerical solutions, and an analytic approximation based on a rescaling of the concentration boundary layer. The discussion of the analytic approximation is self-contained and does not assume any familiarity with the associated boundary layer methods in fluid mechanics. The pedagogical exposition is intended so that an advanced undergraduate student with exposure to boundary value problems common in a quantum mechanics or electricity and magnetism course can follow the derivation.

In what follows we consider the flow profile around a single sphere due to deformations of its surface, which is a common model for active particles felderhof2015stokesian (). This model is chosen for several reasons. First, it provides a model to mimic the flow around a spherical cell swimming through a liquid. Second, since the solution for the flow profile is known, we can make progress on an analytic approximation to the associated advection-diffusion equation. Our analysis of the concentration profile focuses on the first-passage probability of a particle in the extracellular fluid being absorbed by the cell redner2001guide (). The first-passage probability reveals the theoretical maximum efficiency for the capture of small molecules by the cell berg1977physics ().

## Ii Fluid Flow

Consider a spherical cell of radius in a fluid. We assume the fluid is infinite in extent and also at rest at infinity. In the case of present interest, the Reynolds number corresponding to the flow around the cell is generally much less than one. At small Reynolds numbers, the magnitude of the inertial terms in the Navier-Stokes equation become much less than the magnitude of the viscous terms. The disparity in magnitude allows for a linearization of the Navier-Stokes equation into the Stokes equation, which neglects the inertial terms.

The flow velocity satisfies the Stokes equation

(1) |

where is the pressure profile and is the dynamic viscosity of water. Since the speed of the flow will be much less than the speed of sound in water, the fluid is considered to be perfectly incompressible and satisfy the incompressibility condition

(2) |

The fluid is set into motion through axially symmetric deformations of the spherical surface of the cell. At small Reynolds number the fluid motion is then determined by the no-slip boundary condition at the surface. The details for the calculation of the flow velocity can be found in felderhof2015stokesian (). Defining the dimensionless length , in the instantaneous rest frame of the cell, the dimensionless fluid velocity is where

(3) |

and is the translational velocity of the cell. The defining relations for and are

(4) | |||

(5) |

with the Legendre polynomial and the associated Legendre function.

## Iii Numerical Model

With a solution for the flow profile in hand, we now turn to the second problem, determining the concentration profile from the advection-diffusion equation. Our eventual aim is to find the first passage probability of a dissolved molecule released at position () through the inner boundary at . To accomplish this, we solve the advection-diffusion equation using the known velocity profile. In this section we first discuss a numerical solution of the problem. Later we will compare our results to an analytic approximation.

For an incompressible fluid, the advection-diffusion equation reads as

(6) |

where is the diffusion coefficient for the molecule of interest, is the molecular concentration, and is time. We define the dimensionless concentration and dimensionless time . In what follows we will work with a time independent flow velocity by considering a temporal average over a period of the fast swimming motion, in which case the Péclet number can be defined as . The dimensionless form of Eq. (6) in spherical polar coordinates is

(7) |

Building the numerical model is generally simpler and less prone to sporadic oscillations in the solution if a grid with uniform spacing is used. This is achieved by casting Eq. (7) in Cartesian coordinates to obtain a numerical solution. Replacing the derivatives in Eq. (7) with their finite-difference counterparts yields

(8) |

Here is the value of the concentration at the location indexed by at the time indexed by with the time step and the distance between grid points. We work on a grid 10 units wide with 200 equidistant points along each axis. Here a unit corresponds to the dimensionless distance equal to one cell radius. The average velocity field is computed first and stored in memory for each grid point. The fluid model used to obtain our numerical solutions, the squirming swimmer, is defined in section IV. B. of felderhof2015stokesian (). The only non-zero multipole moments are , , and . We use an explicit finite difference scheme where the concentration at time is computed from the concentration at time . This is continued until the percent change in the concentration field is no larger than 0.01% per iteration at any grid point: indicating that the concentration field has reached steady state. The first-passage probability is then determined as the flux through the inner absorbing boundary divided by the sum of the flux through the inner and outer boundaries.

## Iv Analytical Model

We now consider an approximate analytical approach to the problem. The basic idea is to use methods common to boundary layer problems acrivos1960solution (); acrivos1962heat (); acrivos1965asymptotic (); goddard1966asymptotic () to find the leading order solution of Eq. (7) with the boundary conditions given. Readers familiar with this literature will recognize the general strategy, but our exposition is intended to be entirely self-contained so that a physics student with no previous exposure to the fluids literature but some exposure to boundary value problems (at the level common in an undergraduate electricity and magnetism or quantum mechanics course) can follow along. The general structure for the calculation is as follows. A first change of variables () defines the boundary layer equation for our problem. A second change of variables (the similarity transformation ) gives us an effective radial equation. This equation is then solved using methods similar to those students are familiar with from solving boundary value equations (Green’s function).

Defining the small parameter , we first rescale the radial variable to stretch out the boundary layer. The dimensionless time is rescaled as , but we do not rescale the angular variables and . Dominant balance rnmg () determines the value of exponents and . By making this choice of exponents, in a perturbative expansion for the concentration, , the equation governing will contain temporal, advective, and diffusive terms. The physical thickness of the concentration boundary layer is . The fluid model is the same one used to obtain our numerical solutions. Inserting the perturbative expansion into Eq. (7) and collecting terms of the same order in , the result is a system of coupled equations for the . Defining , the lowest order governing equation is

(9) |

A common method in this type of boundary layer problem nuss (); anu () is to define similarity variables and , where encapsulates the angular dependence of the boundary layer. Making this choice, the following relationships are needed for the change of variables: , , , . This transformation isolates the angular dependence and defines the concentration boundary layer equation for the problem,

(10) |

Provided the term in large parenthesis is equal to a constant , the equation is transformed into an effective radial equation. To proceed we define the time-integrated concentration

(11) |

The equation for becomes

(12) |

Before solving this equation for the first-pasage probability, we discuss the solution to the associated angular equation for . In what follows we make the choice , in which case we obtain a first-order differential equation for the variable ,

(13) |

where

(14) | |||||

(15) |

This equation can be solved with an integrating factor. To do so we calculate which gives

(16) |

The desired solution is then

(17) |

where is an integration constant.

We are now in a position to proceed with the calculatin of the first-passage probability. Consider a spatial domain where all molecules released are eventually captured with probability one. The concentration is defined in the spherical shell between two perfectly absorbing surfaces, the first at the surface of the cell , and a second at some prescribed distance . Since all molecules are eventually absorbed, . If a molecule is released in the extracellular fluid at position , the corresponding initial condition is .

Writing the initial condition in terms of , our effective radial equation, Eq. (12), becomes

(18) |

To solve Eq. (18), note that the two independent solutions to the homogeneous equation are a constant and . Using perfectly absorbing boundary conditions at and , the solution for is

(19) |

Here () is the smaller (larger) of and . To determine the constant , we integrate both sides of the governing equation to determine the discontinuity in the first derivative of ,

(20) |

where . The factor of in the integration to determine the discontinuity was incorrectly omitted in urchin (), it is needed to ensure the proper normalization,

(21) |

In terms of the original variables, the first-passage probability is calculated as

(22) |

The result can be written in terms of the boundary layer variable and the time-integrated concentration as

(23) |

Performing the angular integration we arrive at

(24) |

Moving the outer absorber out to infinity and writing the result back in terms of yields the final result,

(25) |

Note that the final result is properly normalized as . The result has certain structural similarities to a much simpler calculation, the first-passage probability in the case of pure diffusion, for which redner2001guide (). In words, the effect of fluid flow is an exponential supression of the first-passage probability as compared with the purely diffusive case.

## V Results

To make a comparison between the analytical approach and the numerical solution for a variety of point source locations, we numerically calculate the average first-passage probability, from Eq. (25). Here, the integration constant in the equation for is zero, . The full set of angular results from the numerical calculation are fit to a function of the form,

(26) |

using a non-linear least squares routine in Python.
Figures 3 and 4 show a comparison of the numerical result for the first-passage probability and the analytic approximation. The result for the first-passage probability is in good agreement with the numerical result, especially at large Pe (see Fig. 3). This is reassuring as the perturbation program is constructed based on a small parameter Pe. As Pe is reduced the agreement becomes less quantitatively accurate (see Fig. 4), but contintues to capture the qualitative behavior, even for source locations far outside of the concentration boundary layer.

To make contact with the motivations for the study discussed in the introduction, note that the cost and effectiveness of the toxin transport system is an active area of experimental research cole2013cost (). A physical microvilli length of m corresponds to a source location of in Figs. 3 and 4. For a cell in a flow field (either generated by its own swimming motion or an environmental flow), designing a smart toxin disposal system would probably only require that the cell have knowledge of the upstream direction. In the biophysical context, the receptors responsible for effluxing toxin molecules to the exerior of the cell do not chemically modify the toxin. As a result, there is a potentailly costly scenario in which effluxed molecules are reabsorbed and have to be discarded again, which is known as futile cycling. The cell could see significant gains in efficiency be either preferentially activating transport receptors downstream or actively transporting molecules tagged for export to downstream receptors. A future experiment designed to monitor receptor activity in a controlled flow environment, perhaps a microfluidic chamber, might be able to uncover whether a smart disposal mechanism similar to this has evolved naturally.

## Vi Conclusion

Motivated by a mass transport problem during embryonic development, we considered an active particle model for a swimming cell, the spherical squirmer. Within the context of the model, we determined the first-passage probability for an advection-diffusion equation at small Reynolds number and large Péclet number. Numerical approaches to solving the advection-diffusion equation based on explicit finite-differencing were compared to an analytic approximation based on a rescaling of the concentration boundary layer. For large Péclet number we find quantitative agreement between the approaches. As the Péclet number is reduced, the analytic approximation continues to capture the qualitative behavior of the first-passage probability, but deviates somewhat from the numerical results. The regime of validity for the analytic approximation might be improved by continuing the perturbation program beyond the zeroth order. This is a potential direction for future research.

## Vii Acknowledgements

This work was supported by faculty startup grant U035843 from the College of Arts, Sciences, and Letters at the University of Michigan-Dearborn.

## References

- (1) B. J. Cole, A. Hamdoun, and D. Epel, “Cost, effectiveness and environmental relevance of multidrug transporters in sea urchin embryos,” Journal of Experimental Biology, vol. 216, no. 20, pp. 3896–3905, 2013.
- (2) M. B. Short, C. A. Solari, S. Ganguly, T. R. Powers, J. O. Kessler, and R. E. Goldstein, “Flows driven by flagella of multicellular organisms enhance long-range molecular transport,” Proceedings of the National Academy of Sciences, vol. 103, no. 22, pp. 8315–8319, 2006.
- (3) K. Whalen, A. M. Reitzel, and A. Hamdoun, “Actin polymerization controls the activation of multidrug efflux at fertilization by translocation and fine-scale positioning of abcb1 on microvilli,” Molecular biology of the cell, vol. 23, no. 18, pp. 3663–3672, 2012.
- (4) X. Song. PhD thesis, Indiana University, 2008.
- (5) S. Vogel, Life in moving fluids: the physical biology of flow. Princeton University Press, 1994.
- (6) T. Gökirmak, J. P. Campanale, L. E. Shipp, G. W. Moy, H. Tao, and A. Hamdoun, “Localization and substrate selectivity of sea urchin multidrug (mdr) efflux transporters,” Journal of Biological Chemistry, vol. 287, no. 52, pp. 43876–43883, 2012.
- (7) N. A. Licata and A. Clark, “Fluid flow enhances the effectiveness of toxin export by aquatic microorganisms: a first-passage perspective on microvilli and the concentration boundary layer,” Physical Review, vol. 91, no. 012709, 2015.
- (8) B. Felderhof, “Stokesian spherical swimmers and active particles,” Physical Review E, vol. 91, no. 4, p. 043018, 2015.
- (9) H. C. Berg and E. M. Purcell, “Physics of chemoreception.,” Biophysical journal, vol. 20, no. 2, p. 193, 1977.
- (10) G. K. Batchelor, An Introduction to Fluid Dynamics. Cambridge University Press, 1967.
- (11) W. S. Janna, Introduction to Fluid Mechanics. PWS Publishing Company, 1993.
- (12) A. Acrivos, “Solution of the laminar boundary layer energy equation at high prandtl numbers,” Physics of Fluids (1958-1988), vol. 3, no. 4, pp. 657–658, 1960.
- (13) A. Acrivos and T. D. Taylor, “Heat and mass transfer from single spheres in stokes flow,” The Physics of Fluids, vol. 5, no. 4, pp. 387–394, 1962.
- (14) A. Acrivos and J. Goddard, “Asymptotic expansions for laminar forced-convection heat and mass transfer,” Journal of Fluid Mechanics, vol. 23, no. 2, pp. 273–291, 1965.
- (15) J. Goddard and A. Acrivos, “Asymptotic expansions for laminar forced-convection heat and mass transfer part 2. boundary-layer flows,” Journal of Fluid Mechanics, vol. 24, no. 2, pp. 339–366, 1966.
- (16) R. H. P. John C. Tannehill, Dale A. Anderson, Computational Fluid Mechanics and Heat Transfer. Taylor and Francis, second ed., 1997.
- (17) O. Zikanov, Essential Computational Fluid Dynamics. John Wiley and Sons, 2010.
- (18) S. Biringen and C.-Y. Chow, Introduction to Computational Fluid Mechanics by Example. John Wiley and Sons, 2011.
- (19) S. Redner, A guide to first-passage processes. Cambridge University Press, 2001.
- (20) L.-Y. Chen, N. Goldenfeld, and Y. Oono, “Renormalization group and singular perturbations: Multiple scales, boundary layers, and reductive perturbation theory,” Phys. Rev. E, vol. 54, pp. 376–394, Jul 1996.
- (21) V. Magar, T. Goto, and T. J. Pedley, “Nutrient uptake by a selfpropelled steady squirmer,” Q. J. Mech. Appl. Math., vol. 56, pp. 65–91, 2003.
- (22) V. Magar and T. J. Pedley, “Average nutrient uptake by a self-propelled unsteady squirmer,” Journal of Fluid Mechanics, vol. 539, pp. 93–112, 9 2005.