We consider the configuration shown in Fig 1, similar to the one addressed in Ref. [23]. Both the cell and the nucleus are modelled as stationary rigid spheres, of radius R_c and R_nuc, respectively. A Cartesian coordinate system is affixed at the centre of the cell, with gravity acting along the negative z-direction. Without loss of generality, the centre of the nucleus is situated in the y = 0 plane at a distance of eR_c from the cell centre, at an angle θ_e from the x-axis. From geometry, it is clear that one must have eR_c + R_nuc ≤ R_c. The quantity eR_c cos θ_e = e_xR_c denotes the eccentricity, along the x-axis, of the nucleus centre relative to the cell centre. Similarly, eR_c sin θ_e = e_zR_c denotes the nucleus eccentricity along the z-axis.

We assume that the cellular membrane is maintained at a constant temperature T_mem, whereas the surface of the nucleus is maintained at a constant temperature T_nuc. Following experimental evidence [5, 9], we take the nucleus to be warmer than the cell membrane (i.e. T_nuc > T_mem). Note that while we are modelling flows due to temperature gradients between the nucleus and the cell membrane, a change in the size (R_nuc) and position (e, θ_e) of the inner sphere can also provide an approximation of convection caused by other warm organelles inside the cell.

With cells that are typically tens of microns in size [30, 31], intracellular fluid flow is expected to be dominated by viscous forces and thus governed by classical microhydrodynamics [32]. Following Ref. [23], we model the aqueous cytoplasm as a Newtonian viscous liquid with physical properties similar to water at a reference membrane temperature, T_mem [30, 31].

We further approximate the viscosity of the cytoplasm as uniform; although viscosity is known to vary with temperature, including this effect is known to only give rise to a small correction to the flow [33]. The driving of the flow is thus solely due to the density of the cytoplasm, which varies spatially in response to thermal gradients and thus necessitates solving for the temperature distribution within the cell. The temperature is governed by an advection-diffusion equation, denoting a balance between heat transfer by fluid flow and by diffusion (conduction). In the next section, we formalise these statements and present the equations governing the flow and temperature fields inside the cell.

In the absence of inertia, the flow field, u′(x′) at a position x′ inside the cell, is governed by the incompressibility constraint (or the continuity equation), ∇ ′ · u ′ = 0 , (1) and the Stokes equations, - ∇ ′ p ′ + η ∇ ′ 2 u ′ + ρ ( T ′ ) g = 0 , (2) where the primes signify dimensional variables.

In the above equation, p′ is the pressure in the fluid, η is the constant dynamic viscosity, g = −gi_z is the acceleration due to gravity (i_z is the unit vector in the z-direction), and T′ is the temperature in the cell. We denote the temperature-dependent density of the fluid by ρ(T′). For sufficiently small temperature differences T′ − T_mem, we model this temperature-dependence with a standard linear relationship, ρ ( T ′ ) = ρ 0 [ 1 + β ( T mem - T ′ ) ] , (3) where β is the thermal expansion coefficient, with units of K⁻¹. It is important to note here that we are employing the incompressibility constraint, Eq (1), yet allowing the density to vary with temperature. This is a classical (and validated) method in hydrodynamics called the Boussinesq approximation [33]. For modest changes in temperature, we can neglect the temperature dependence of material properties everywhere except in the gravitational force term in Eq (2), when this buoyancy term is the driving mechanism for the flow (as in the present problem). This provides a first approximation to the fluid flow, which should be quantitatively accurate for temperature differences less than ∼10 K–20 K [33]. Mathematically, this may be formalised via a perturbation expansion of the full governing equations in the limit of small temperature changes [33].

The governing equations, Eqs (1) and (2), must be complemented by appropriate boundary conditions for the velocity. Since we are modelling the nucleus and the cell as stationary, rigid spheres, the fluid velocity must vanish at both these surfaces, i.e. u ′| nuc = 0 , u ′| mem = 0 , (4) where the sub-scripts ‘nuc’ and ‘mem’ denote the cell nucleus and membrane, respectively. It is clear from Eq (2) that solving for the fluid flow requires knowledge of the temperature. The steady-state temperature field, T′(x′), is governed by an advection-diffusion equation, u ′ · ∇ ′ T ′ = α ∇ ′ 2 T ′ , (5) subject to the boundary conditions, T ′| nuc = T nuc , T ′| mem = T mem . (6) In Eq (5), α is the constant thermal diffusivity of the cytoplasm (assumed to be water in the present study).

It is important to note that for the geometry described in Fig 1, Eqs (1) to (6) cannot admit a quiescent solution (u′ = 0) as long as there exists a temperature contrast between the two surfaces. For any T_nuc ≠ T_mem, the temperature distribution T′ (x′) will be non-uniform and the resulting temperature gradients (and hence the density gradients) will not be aligned with gravity inside the cell. Thus, no hydrostatic pressure distribution can be found that balances the gravitational forcing everywhere, and the system cannot stay in equilibrium without fluid flow. This is fundamentally different from the classical Rayleigh–Bénard convection, where gravity is parallel to an imposed temperature gradient and fluid flow emerges as a result of an instability above a critical imposed temperature gradient [26].

We next render the equations dimensionless; using the characteristic temperature difference ΔT = T_nuc − T_mem, we may define the reference scales for length, velocity and pressure, respectively, as ℓ ref = R c , u ref = ρ 0 β Δ T g ℓ ref 2 η , p ref = η u ref ℓ ref . (7)

We also define a re-scaled temperature, Θ = T ′ - T mem T nuc - T mem . (8)

Using Eqs (7) and (8) in the equation for the flow, Eq (2), yields now in dimensionless form (without primes) - ∇ P + ∇ 2 u + Θ i z = 0 , ∇ · u = 0 , (9) subject to the dimensionless boundary conditions u| nuc = 0 , u| mem = 0 . (10) The pressure P in Eq (9) is now a modified dynamic pressure, given in dimensionless form by P = p + (βΔT)⁻¹ i_z ⋅ x. Similarly, we can re-write the thermal problem, Eq (5), in terms of the normalised temperature Θ(x) as u · ∇ Θ = 1 Pe t ∇ 2 Θ , (11) with the normalised boundary conditions, Θ| nuc = 1 , Θ| mem = 0 . (12) Importantly, in Eq (11), Pe_t ≡ u_refR_c/α is the (dimensionless) thermal Péclet number, a ratio of the rate of heat transfer by fluid flow to the rate of heat transfer by diffusion.

The typical values of the parameters used for calculating the reference scales in our problem are given in Table 1, where, just like Ref. [23], we have assumed that the cytoplasm shares the material properties of water at temperature T_mem = 310 K (≈ 37°C). This is a reasonable assumption for both the physical [30, 34] and thermal [35, 36] properties of the cytoplasm. We note, however, that in some instances the effective viscosity of the cytoplasm can be higher than that of water, due to a high concentration of macromolecules [30]; we will discuss the implications of this disparity in Sec. 4.2 and argue that it does not alter our central result. The size (radius) of the cell in our study corresponds roughly to the typical COS-7 [9] or HeLa [10] cells used in intracellular thermometry studies, while we consider a wide range of sizes and positions of the nucleus, to reflect the diversity in the nucleus’s placement inside biological cells, depending on the cell type and the stage of the cell cycle [27, 28]. Here too, we note that the main conclusions of this study remain independent of the details of nuclear size and position inside the cell (see Appendix B).

A quantitative description of natural convection in the model cell requires us to solve Eqs (9)–(12) for Θ and u. For a prescribed thermal Péclet number (set by the material properties of the cytoplasm), the solution depends only on the size and position of the nucleus. These quantities are captured in three dimensionless numbers: κ = R_nuc/R_c, the ratio of the nucleus’s radius to that of the cell, and (e_x, e_z), the eccentricity of the nucleus centre along the x- and z-axes.

Depending on the specific values of (e_x, e_z) and Pe_t, one can solve Eqs (9)–(12) using full numerical simulations, semi-analytical methods, or even analytically to obtain an explicit expression for the temperature and flow, (Θ, u). In following sections, we describe solutions obtained using each of these methods to provide rigorous estimates of the flow inside the cell.

We first present dimensional results of the temperature and flow fields inside the cell, to gain intuition about the strength of the temperature-gradient-driven flow. We show in the top row of Fig 2 one specific solution, for a radius ratio κ = 0.43, eccentricity values e_x = e_z = 0.25, and a nucleus-to-membrane temperature difference of ΔT = 1 K, evaluated numerically using the finite-elements-based software COMSOL; this is the benchmark geometry used in Ref. [23]. The values of the other physical parameters are stated in Table 1 and are the same as in Ref. [23], allowing us to compare our results directly. The results of Ref. [23] are reproduced in the bottom row of Fig 2. For the flow field, we plot the horizontal (i.e. perpendicular to gravity) velocity u x ′ and the vertical (i.e. upward) velocity u z ′ at the y = 0 mid-plane of the cell. Since this is a plane of symmetry, the velocity component normal to it, i.e. u y ′, is identically zero.

An important implication of analysis in axisymmetric geometry is the symmetry of the flow field around the axis between the nucleus and the centre of the cell. We may exploit this symmetry to solve Eqs (9) and (13) semi-analytically in a bi-spherical coordinate system, (ξ, χ, ϕ), shown in Fig 4(a), where the nuclear membrane and the cell membrane can be represented by distinct coordinate surfaces [32]. Axisymmetry allows the flow field to be represented in terms of a single Stokes streamfunction, Ψ^b(ξ, χ), which can be expressed as a linear superposition of harmonic functions [32]. The technical aspects of the solution are detailed in Appendix D. A sample solution in the axisymmetric case is shown in Fig 4(b), where the z-axis is the axis of symmetry and the geometric parameters are κ = 0.43, e_x = 0, e_z = 0.3 (left: temperature; right: vertical velocity). We can also solve the axisymmetric problem using finite-element COMSOL simulations. In Fig 5, we compare the results obtained from the two different methods, by plotting the dimensionless maximum (max. {u_z}) and minimum (min. {u_z}) vertical velocities inside the cell, as a function of the eccentricity e_z. We provide further validation of the temperature and velocity fields in Fig 16 in Appendix D. The excellent matching between the two sets of results validates the non-axisymmetric simulations from Sec. 3.

We see that the results are symmetric about e_z = 0, which is a result of the structure of Eqs (9) and (13). The temperature distribution in the cell is reflected about the z = 0 plane as e_z goes from negative to positive (see Fig 6(a) and 6(c)). Since the fluid flow is directly proportional to the local temperature difference, regions in the cell corresponding to the same thermal environment in the two geometries also display the same vertical velocity (Fig 6(b) and 6(d)). As a result, the location of the strongest upwelling flow with respect to the nucleus changes: from being above the centre of the nucleus for e_z < 0, to shifting below the nucleus centre for e_z > 0.

In summary, assuming an axisymmetric geometry allowed us to solve for the fluid flow via a semi-analytical approach that is less computationally intensive than the finite-element simulations. By performing two independent sets of calculations in the axisymmetric limit (finite-element simulations and the Stokes streamfunction analysis of this section) we have confirmed the velocity magnitudes that were predicted numerically in Sec. 3. In the next section, we perform a third analysis, in which we consider the limit where the nucleus is located at the centre of the model cell (i.e. e_x = e_z = 0). This particular geometry enables us to obtain an exact, analytical solution for the temperature distribution and the flow field.

When the nucleus is concentric with the cell, one can solve the problem fully analytically in a spherical coordinate system, (r, θ, ϕ) [25]. Since we still have an axisymmetric setup, we can once again represent the flow field in terms of a streamfunction, Ψ^s(r, θ), where r is the radial separation of any point, measured from the centre of the cell (and the nucleus), and θ is the polar angle measured from the positive z-axis [32]. We discuss in the main text the final results, i.e. the analytical expressions for the temperature and the fluid velocity; all details of the solution methodology are provided in Appendix E. The temperature is found to be radially isotropic (i.e. it only depends on r) and is given by Θ ( r ) = κ 1 - κ ( 1 r - 1 ) . (14)

As for the flow field, we derive in Appendix E a very simple representation for the streamfunction Ψ^s(r, θ) as Ψ s ( r , θ ) = f 1 ( r ) sin 2θ , (15) where f 1 ( r ) = r 3 8 κ 1 - κ + c 1 r + c 2 r + c 3 r 2 + c 4 r 4 , (16) with the constants, c_i, being functions solely of the radius ratio κ given by c 1 = κ 2 ( 3 κ 2 + 4 κ + 3 ) 8 ( 1 - κ ) ( 4 κ 2 + 7 κ + 4 ) , c 2 = - κ 4 8 ( 1 - κ ) ( 4 κ 2 + 7 κ + 4 ) , c 3 = - κ ( κ + 1 ) ( κ 2 + 3 κ + 1 ) 4 ( 1 - κ ) ( 4 κ 2 + 7 κ + 4 ) , c 4 = - κ ( κ + 1 ) 4 ( 1 - κ ) ( 4 κ 2 + 7 κ + 4 ) . (17)

Once the streamfunction is known, the fluid velocity u may be obtained explicitly as u = ( u r , u θ ) = ( - 1 r 2sinθ ∂ Ψ s ∂ θ , 1 rsinθ ∂ Ψ s ∂ r ) , (18) where u_r is the radial component and u_θ the tangential component of the fluid velocity. Substituting Eqs (15) in (18) yields, u r ( r , θ ) = - 2 f 1 ( r ) r 2 cosθ , u θ ( r , θ ) = 1 r d f 1 d r sinθ . (19) We thus have a closed-form, exact solution to the temperature-gradient-driven natural convection problem in the concentric case, for arbitrary radius ratios κ = R_nuc/R_c. A simple transformation from spherical to Cartesian coordinates then yields the vertical velocity u_z and the horizontal velocity u_x (see Eq (58) in Appendix E). In Fig 7, we plot these velocity components (left column) and compare them with the corresponding COMSOL simulation results (right column). We can see an excellent agreement between the exact solution and our computations; quantitatively, the average relative error between the two solutions is less than 2%, which provides another validation of the simulations of Sec. 3 and confirms that the velocity predictions therein are accurate.

The results in Eqs (15) to (19) give us an exact representation of the flow inside the cell, in this simple geometrical case. How relevant are these predictions for arbitrary locations of the nucleus? Further computations shown in Fig 13 (Appendix B) demonstrate that the strength of the flow depends only weakly on the position of the nucleus. Therefore, the exact analysis of this section captures all essential aspects of this hydrodynamics problem.

We have now validated our numerical simulations using two independent analyses, which gives us confidence in our numerical results from Sec. 3. We thus conclude that the flow strengths that we report are a more accurate estimate of intracellular buoyancy-driven flows than the results of Ref. [23].

We model the solute as a passive scalar that diffuses and is advected by the flow calculated in Sec. 3. At steady state, the dimensionless concentration of the scalar, C, is governed by the advection-diffusion equation u · ∇ C = 1 Pe s ∇ 2 C , (21) where Pe_s = u_refR_c/D is the Péclet number of the solute (the ratio of its advective to its diffusive transport), where D is the solute’s molecular diffusivity. For the sake of simplicity, we assume that the scalar is produced at the nucleus and absorbed at the cell membrane, so we have the boundary conditions C nuc = 1 , C mem = 0 . (22) Typical molecular diffusivities of cellular matter can be as low as 0.01 μm² s⁻¹ [30, 38], which corresponds to Pe_s ≈ 500. Based on this preliminary scaling approach, one may then expect the cell-scale natural convection described previously to significantly affect solute transport.

To investigate the extent to which this is true, we solve Eqs (21) and (22) for the configuration and parameters of Sec. 3 (κ = 0.43, e_x = e_z = 0.25, ΔT = 1K) using finite-element COMSOL simulations. The resulting solute distributions on the y = 0 mid-plane are shown in Fig 9. Surprisingly, the effect of natural convection is so weak that even for Pe_s as large as 500 (Fig 9(a)), the solute distribution is almost the same as that when the solute transport is purely diffusive, i.e. when Pe_s ≪ 1 (Fig 9(c)). One can also visualise the time evolution of the concentration, starting from a uniform initial condition C(t′ = 0) ≡ 0 and then suddenly changing the concentration at the nucleus surface to C_nuc = 1 at t′ = 0.1 s; this is done in Appendix C. As shown in Fig 14 there, the concentration distribution becomes nearly independent of Pe_s when the time is normalised by the diffusive time-scale t d ′ = R c 2 / D, which is a classic signature of isotropic diffusion. This further confirms that the transport mechanism inside the cell is mostly diffusive.

While the maximum local Péclet number Pe max ℓ captures a local estimate of mass transfer enhancement by cytoplasmic convection, a fairer quantification of this enhancement at the whole-cell level involves calculating the total absorption flux across the nuclear membrane, and comparing it to the limit when the solute emitted by the nucleus instead undergoes pure diffusion. If S denotes the surface of the nucleus, then the dimensionless diffusive flux through the surface is given by the integral ∫_S −n ⋅ ∇C dS, while the dimensionless advective flux, ∫_S n ⋅ u C dS, vanishes due to impermeability of the nuclear surface (see Eq (4)). The influence of bulk solute advection on the mass transfer rate across the nucleus is quantified by defining a normalised flux, J = ∫ Sn · ∇ C d S ( ∫ Sn · ∇ C d S )| Pe s ≡ 0 , (23) where the denominator is the solute flux through the nucleus when the solute transport is purely diffusive (i.e. when Pe_s ≡ 0 in Eq (21)). In the mass transfer literature, the quantity J is usually referred to as the Sherwood number [33] and is expected to be a function of the cell geometry and the solutal Péclet number, Pe_s.

If bulk advection of solute were to significantly impact the solute transfer from the nucleus to the cell membrane, then we would expect J to be greater than unity. We plot in Fig 10(c) the dependence of this normalised flux on the diffusivity of the solute (D) for various temperature differences (ΔT). The solute diffusivity D spans four orders of magnitude, covering cellular material ranging from large vesicles (D ∼ 10⁻² μm² s⁻¹ [39, 40]), to proteins synthesised in the ribosomes (D ∼ 1 μm² s⁻¹ to 10 μm² s⁻¹ [30, 41]), to metabolites such as ATP, ADP and calcium ions (D ∼ 100 μm² s⁻¹ [42, 43]).

For the lowest ΔT value considered, the convective flows do not provide any advantage in driving the solute. Even for ΔT as large as 10 K (a large value not supported by experimental data) and solute diffusivities as low as 0.01 μm² s⁻¹ (relevant for the largest vesicles in the cell), the solute removal from the nucleus is only amplified by about 15% compared with the case of pure diffusion.

As an example, consider the protein insulin, whose diffusivity in living cells is around 1 μm² s⁻¹ [30]. The result in Fig 10(c) then tells us that J ∼ 1 for any range of experimentally observed intracellular temperature differences. Most cellular proteins and macromolecules have diffusivities that are too large (ranging from 1 μm² s⁻¹ to 10 μm² s⁻¹) to be overcome by natural convection. Hence, their motion through the cell likely remains diffusive. Vesicles in cells can have very low diffusivities, around 10⁻² μm² s⁻¹ [39, 40] but, despite this, for the commonly measured temperature differences within cells, their transport enhancement due to convective flows is expected to remain, at best, very modest, and most likely, completely negligible.

In this work, we re-considered the problem of intracellular natural convection potentially induced by a temperature difference ΔT between the cell’s nucleus and membrane. Specifically, we used numerical simulations to show that intracellular flows are expected to be on the order of ∼ 10⁻³ μm s⁻¹ in magnitude for a nucleus-to-membrane temperature difference of 1 K. In the process, we discovered an important disparity from previous results that had predicted ten-fold stronger flows for the same geometry and temperature differential [23].

We explored further this discrepancy, and confirmed it, by performing independent calculations of the flow strength in an axisymmetric geometry (Sec. 4) and in a concentric geometry (Sec. 5). Importantly, in the latter case, we provided an exact expression for the flow inside the cell. All our analyses converged to a similar value for the flow strength, which was indeed ∼10 times weaker than that computed in Ref. [23]. These results highlight the importance of detailed flow calculations, since simple scaling analysis alone suggests fluid velocities almost two orders of magnitude larger (u_ref ∼ 0.5 μm s⁻¹), which would lead one to believe, incorrectly, that these convective flows are very strong.

The transport of materials within a cell is essential to its normal function, from vesicle transport for structural upkeep [44], signal transmission via proteins [41] and organelle transport during cell division [45, 46] to inter-organelle interactions [47] and the maintenance of nutrient/metabolite gradients [48]. The movement of cellular matter is accomplished by a variety of mechanisms, including molecular diffusion, active transport by motor proteins and advection due to cytoplasmic flows [37]. A second aim of our work was to investigate whether temperature-gradient-driven flows could contribute usefully to advective material transport inside the cell. This would be particularly relevant for the motion of cell constituents with low diffusivities (large organelles and vesicles), which necessarily require alternative mechanisms to traverse the cell. Towards this, we numerically simulated the mass transfer of a chemical species by cytoplasmic convection, and showed that convection does not contribute significantly to the transport of material within a cell beyond what is achieved by pure diffusion. We quantified the advection-induced enhancement in averaged mass-flux of a chemical species released from the nucleus and absorbed at the cell membrane, as a function of the prescribed temperature difference and the diffusivity of the species. Only in situations that are biologically unrealistic (temperature differences of 10 K and very small diffusivities of 0.01 μm² s⁻¹, associated with cellular vesicles) could the averaged mass-flux see a modest increase above that obtained with purely diffusive mass transfer.

Thus, while a cell’s thermal environment is important for its metabolism [3, 49] and survival [50], the flows generated by temperature gradients, in most realistic scenarios, have negligible effect in improving intracellular mass transfer. It seems therefore that these flows do not noticeably impact important cell processes like protein delivery, signal propagation, and organelle and metabolite transport.

Once again, this conclusion becomes apparent only after a detailed calculation of solute transport in the bulk, as done in Sec. 6. A simple scaling analysis yields characteristic solutal Péclet numbers Pe_s = u_refR_c/D ≈ 500 that are much larger than unity, which, in the absence of detailed numerical calculations, would erroneously suggest a significant natural-convection-induced advective contribution to species transport.

One of the reasons for the aforementioned flows being so weak is that we have modelled the nucleus surface and the cell membrane as rigid and non-slipping, i.e. the fluid’s tangential velocity vanishes at both these surfaces. As a result, there are significant viscous stresses in the domain that resist fluid motion. A different model could allow tangential motion along the cell surface due to slip or mobility of the membrane.

Our model simplifies the cytoplasm and assumes that it is effectively a homogeneous continuum with Newtonian properties. It is obviously known to be more complex, consisting of a polymeric and dynamic cytoskeleton embedded in a viscoelastic, gel-like fluid that flows past freely suspended cell organelles [31].

Our analysis also considers a simplified description of thermal diffusion inside the cell. Inherent in Eq (5) is the assumption that the thermal diffusivity of the cytoplasm is isotropic, i.e. that heat diffuses at a constant rate along all directions. Since the cytoplasm is heterogeneous, the thermal diffusivity within it is expected to vary spatially [51].

Our calculations show that cytoplasmic flows due to naturally occurring temperature gradients within a cell are not very effective in driving cellular material transport. However, artificially induced thermal convection has been recently proposed as a means to accelerate cell assembly for biomedical assays [52]. Light-absorbing particles are arrested in optical traps, are heated in the process and thus induce thermal flows that drive cell accumulation around the particles. One could envision the generation of such flows within a cell as well, for example, by laser-assisted heating of metallic nano-particles [53–55]. Our theoretical analysis could thus be adapted to quantify such artificially induced cell-scale convection, to inform future hydrodynamics-based strategies for intracellular object manipulation.

In the main text, we discussed the nature of fluid flow in the vertical plane through the centres of both the cell and the nucleus (the x′-z′ plane). Here we discuss the importance of the other component of the fluid velocity, i.e. the velocity u y ′, in comparison with u x ′ and u z ′. In Fig 11 we show this horizontal velocity u y ′ on two vertical planes x′ = const.—through the centre of the cell and through the centre of the nucleus—along with the vertical velocity u z ′. Similarly, Fig 12 shows u y ′ on five different horizontal planes through the cell, at different heights z′. The cell geometry and the physical parameters in both these figures are the same as in Fig 3. Note that since only the full three-dimensional flow is incompressible, the in-plane velocity in Figs 11 and 12 is not incompressible; this is most apparent in Fig 12, where we can see local two-dimensional sources and sinks of flow. Clearly, the typical magnitude of u z ′ is larger than u y ′ in both Figs 11 and 12 (see also Fig 3). This means that the strongest flows in the cell are vertical (along or opposite to gravity) and justifies our focus on u z ′ in the main text.

In this section, we use numerical simulations to briefly summarise the dependence of the flow on the position of the nucleus inside the cell. The angle θ_e = tan⁻¹(e_z/e_x) shown in Fig 1 quantifies how axisymmetric the system is, with θ_e = π/2 rad (= 90°) representing the axisymmetric case. Results plotted in Fig 13 show that, for all values of the eccentricity e, the flow strength increases monotonically as the extent of axisymmetry decreases. Intuitively, for θ_e = 0, the majority of the fluid is least confined in the direction normal to gravity. Hence, for any non-zero eccentricity/offset, the fluid experiences the least viscous resistance when θ_e = 0. As a result, it attains the largest velocities for a prescribed driving temperature difference. This can be seen in Fig 13(c), where the fluid to the left of the nucleus is heated and attains large upwelling speeds.

The geometry of the bi-spherical coordinate system is shown in Fig 15. The Cartesian coordinates (x, y, z) are expressed in terms of the bi-spherical coordinates (ξ, χ, ϕ) as x = a 1 - χ 2 coshξ - χ cosϕ , y = a 1 - χ 2 coshξ - χ sinϕ , z = asinhξ coshξ - χ , (24) where a = |sinh(ξ_mem)|. The bi-spherical coordinates enable us to represent the nuclear and cell membranes as surfaces with a constant value of the ξ coordinate. These values at the nucleus (ξ_nuc) and the membrane surface (ξ_mem) are functions of the eccentricity e_z and the radius ratio κ, ξ mem = - e z| e z | cosh - 1 ( 1 - κ + e z 2 2| e z | ) , ξ nuc = - e z| e z | cosh - 1 ( 1 - κ - e z 2 2 κ | e z | ) . (25) Note that |ξ_nuc| is always greater than |ξ_mem|, and that both ξ_nuc and ξ_mem are positive for e_z < 0 and negative for e_z > 0. With this geometrical setup, we may now describe the solution of Eqs (9) and (13) in bi-spherical coordinates.

The flow field is axisymmetric and is expressed in terms of the Stokes flow streamfunction, Ψ^b(ξ, χ), as u = u ξ i ξ + u χ i χ , = ( coshξ - χ ) 2 a 2 ( ∂ Ψ b ∂ χ i ξ - 1 1 - χ 2 ∂ Ψ b ∂ ξ i χ ) , (30) where the super-script ‘b’ denotes the coordinates (bi-spherical) in which the streamfunction is defined. The velocity components (u_ξ, u_χ) are directed normal to the (ξ, χ) iso-surfaces (see Fig 15) and are positive in the direction of increasing (ξ, χ). The unit vectors (i_ξ, i_χ) are expressed in terms of the Cartesian coordinate unit vectors (i_x, i_y, i_z) as i ξ = 1 - χcoshξ coshξ - χ i z - 1 - χ 2sinhξ coshξ - χ ( i xcosϕ + i ysinϕ ) , i χ = 1 - χ 2 sinhξ coshξ - χ i z + 1 - χcoshξ coshξ - χ ( i xcosϕ + i ysinϕ ) . (31) The equation governing the streamfunction can be derived from the Stokes Eq (9) as [32, 33] - coshξ - χ a 1 - χ 2 E 2 [ E 2 ( Ψ b ) ] ( - i xsinϕ + i ycosϕ ) = ∇ × ( Θ i z ) , (32) where E² is the differential operator E 2 ( Ψ b ) ≡ coshξ - χ a 2 [ ∂ ∂ ξ ( ( coshξ - χ ) ∂ Ψ b ∂ ξ ) + ( 1 - χ 2 ) ∂ ∂ χ ( ( coshξ - χ ) ∂ Ψ b ∂ χ ) ] . (33)

The general solution to Eq (32) requires us to write Ψ^b as Ψ b ( ξ , χ ) = ( coshξ - χ) - 3 / 2 ∑ n = 0 ∞ ( 1 - χ 2 ) d L n d χ U n ( ξ ) , (34) where U_n(ξ) are unknown functions (velocity modes), that will depend linearly on the temperature modes Θ_n (Eq (28)). If the summations in Eqs (26) and (34) are truncated after N terms, then we need to solve for N functions U_n(ξ), 1 ≤ n ≤ N. Substituting Eqs (34) and (26) into Eq (32) yields ∑ n = 0 ∞ - Γ 3 2 ( n + 1 ) a 4 L n + 1 ( χ ) - χ L n ( χ ) 1 - χ 2 E n b ( ξ ) = ∑ p = 0 ∞ ( 1 - χcoshξ ) d L p d χ Θ p ( ξ ) - ∑ p = 0 ∞ L p ( χ ) 2 ( Θ p ( ξ )coshξ + 2 d Θ p d ξsinhξ ) , (35) where Γ ≡ (cosh ξ − χ), and E n b ( ξ ) is the function E n b ( ξ ) = 1 2 d 4 U n d ξ 4 - ( n 2 + n + 5 4 ) d 2 U n d ξ 2 + ( n 4 2 + n 3 - n 2 4 - 3 n 4 + 9 32 ) U n ( ξ ) . (36)

We next multiply Eq (35) by (1 − χ²)(dL_i/dχ) and integrate over the limits χ = −1 to χ = 1. When this projection is carried out for 1 ≤ i ≤ N, it yields a system of N coupled, linear ordinary differential equations for the velocity functions {U₁(ξ), U₂(ξ), …, U_N(ξ)}. These functions are expressible in terms of the known modal distribution of the temperature Θ_p(ξ). The projection onto (1 − χ²)dL_i/dχ of the right-hand-side of Eq (35) can be simplified using the properties of Legendre polynomials as 2 i ( i + 1 ) 2 i + 1 [ Θ i ( ξ ) - coshξ 2 { Θ i + 1 ( ξ ) + Θ i - 1 ( ξ ) } - sinhξ { d Θ i - 1 / d ξ 2 i - 1 - d Θ i + 1 / d ξ 2 i + 3 } ] . (37)

The projection onto (1 − χ²)dL_i/dχ of the left-hand-side of Eq (35) yields - 2 a 4 ∑ n = 0 ∞ ∫ - 1 1 ( coshξ - χ) 3 ( n + 1 ) { L n + 1 ( χ ) - χ L n ( χ ) } d L i d χ E n b ( ξ ) d χ , (38) which can be written in short-hand as - 2 a 4 ∑ n = 0 ∞ I i n ( ξ ) E n b ( ξ ) , (39) where I i n ( ξ ) = ∫ - 1 1 ( coshξ - χ) 3 ( n + 1 ) { L n + 1 ( χ ) - χ L n ( χ ) } d L i d χ d χ = ∫ - 1 1 i ( n + 1 ) ( coshξ - χ) 3 { L n + 1 ( χ ) - χ L n ( χ ) } χ L i ( χ ) - L i - 1 ( χ ) χ 2 - 1 d χ . (40)

The second line of Eq (40) is obtained by making use of the following recurrence relation for Legendre polynomials, 1 - χ 2 i d L i d χ = L i - 1 ( χ ) - χ L i ( χ ) .

In left-hand-side of Eq (40), the index ‘i’ in I i n denotes the term onto which the projection is made: (1 − χ²)dL_i/dχ, and the index ‘n’ denotes the contribution to the projection from the n^th term in the expansion on the left-hand-side of Eq (35).

The next step is to write out (39) for a discrete number of points, say M, along ξ. Thus, each mode U_n(ξ) is evaluated on M points: {U_n(ξ₁), U_n(ξ₂), …, U_n(ξ_M)}, where ξ₁ = ξ_nuc and ξ_M = ξ_mem. The discrete version of Eq (39) at ξ = ξ_l(1 ≤ l ≤ M) reads - 2 a 4 ∑ n = 0 ∞ I l i n E n b ( ξ l ) , (41) where I l i n is a three-dimensional matrix resulting from the evaluation of I i n ( ξ ) at ξ = ξ_l, i.e. I l i n = I i n ( ξ = ξ l ). If N is the number of modes (or the upper limit of the summation) at which the temperature (Eq (26)) and the streamfunction (Eq (34)) expansions are truncated, then the size of I l i n is M × N × N. At this level of modal resolution/truncation, Eq (35) written for ξ = ξ_l and projected onto ( 1 - χ 2 ) d L i d χ is thus given by - 2 a 4 ∑ n = 0 N I l i n E n b ( ξ l ) = 2 i ( i + 1 ) 2 i + 1 [ Θ i ( ξ l ) - coshξ l 2 { Θ i + 1 ( ξ l ) + Θ i - 1 ( ξ l ) } ] - 2 i ( i + 1 ) 2 i + 1 [ sinhξ l ( d Θ i - 1 / d ξ 2 i - 1 - d Θ i + 1 / d ξ 2 i + 3 )| ξ = ξ l ] . (42)

While the right-hand-side of Eq (42) is explicitly known (see Eq (28)), the derivatives of the functions U_n(ξ) (inherent in definition of E n b ( ξ ); see Eq (36)) in the left-hand-side need to be discretised. In the present work, we use second-order accurate finite differences for this discretisation. This leads to a total of N × M unknowns, U_n(ξ_l), where 1 ≤ n ≤ N and 1 ≤ l ≤ M.

Eq (42) is a fourth-order, linear, coupled ordinary differential equation for the functions U_n(ξ). It is supplemented by four boundary conditions involving velocity components (u_ξ, u_χ) at the inner and outer sphere surfaces (i.e. at the nucleus and the cell membrane), u ξ ( ξ = ξ nuc ) = u ξ ( ξ = ξ mem ) = 0 , u χ ( ξ = ξ nuc ) = u χ ( ξ = ξ mem ) = 0 . (43)

Using Eq (30), one can show that these boundary conditions yield the following equations for U_n(ξ) and dU_n/dξ, U n ( ξ nuc ) = U n ( ξ mem ) = 0 , d U n d ξ| ξ = ξ nuc = d U n d ξ| ξ = ξ mem = 0 . (44)

Eqs (42) and (44) provide the N × M linear equations required to obtain the unknown functions U_n(ξ_l). Inverting this system of equations solves the hydrodynamic problem for the flow field u(ξ, χ) (obtained via Eqs (30) and (34)). A transformation from u(ξ, χ) to u(x, y, z) (using Eq (31)) allows us to plot and compare the velocities obtained using the above methodology with the finite-element COMSOL simulations. This comparison is shown in Fig 16 where we obtain essentially an identical match between the two solutions. Note that in bi-spherical coordinates, we have solved the problem in a grid where the z-coordinate of the cell centre is not given by z = 0 (see Fig 15), but while plotting our results to compare with those from COMSOL simulations, we have shifted the domain along the z-axis such that the cell centre always lies at (x = 0, z = 0).