2. General solutions

We begin by expressing the conservation equations in the two-phase model, which includes momentum conservation equations for the fluid and the network phases as well as the mass conservation of the mixture. The network and the fluid mechanics are coupled through the friction forces that arise from the relative motion of these two phases. Due to the microscopic dimensions, inertial forces can be neglected in determining the intracellular mechanics. With these assumptions, the conservation equations are as follows.

(2.1a)$$\begin{gather} \text{Mass conservation:} \quad \boldsymbol{\nabla}\boldsymbol{\cdot} [\phi \boldsymbol{v}_n+(1-\phi)\boldsymbol{v}_f]=0. \end{gather}$$

(2.1b)$$\begin{gather}\text{Momentum equation for the network:} \quad \boldsymbol{\nabla}\boldsymbol{\cdot} [\boldsymbol{\sigma}_n-\phi p \boldsymbol{I}]-\xi (\boldsymbol{v}_f-\boldsymbol{v}_n)=\boldsymbol{0}. \end{gather}$$

(2.1c)$$\begin{gather}\text{Momentum equation for the fluid:} \quad \boldsymbol{\nabla}\boldsymbol{\cdot} [\boldsymbol{\sigma}_f-\left(1-\phi\right) p \boldsymbol{I}]+\xi(\boldsymbol{v}_f-\boldsymbol{v}_n)=\boldsymbol{0}. \end{gather}$$

Here, subscripts $n$ and $f$ denote network and fluid phases, $\phi$ is the volume fraction of the network phase, $\boldsymbol {\sigma }_{n}$ and $\boldsymbol {\sigma }_{f}$ are the network and fluid stresses, $p$ is the pressure that ensures the incompressibility constraint of the fluid phase, $\xi$ is the friction constant per unit volume that couples the network and fluid dynamics, and $\boldsymbol {I}$ is the identity matrix. The governing equations form a well-posed linear set of PDEs, when the network stress is modelled using linear VE constitutive equations (CEs). We take the volume fraction of the network phase, $\phi$, to remain constant in space and time. Thus the compressibility of the network leads to non-zero divergence of the fluid velocity field: $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v}_f =(\phi /(\phi -1))\,\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v}_n \neq 0$. However, this apparent fluid compressibility is not associated with changes in fluid density, and arises from fluid sources and sinks in response to changes in the volume of the network phase. As such, the Newtonian fluid stress does not include the isotropic term that scales with $\eta _b (\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {v}_f)$, where $\eta _b$ is the bulk viscosity.

We use a general linear VE isotropic CE to describe the traceless component of the fluid stress:

(2.2)\begin{equation} {\boldsymbol{\sigma}}_f = \int_{0}^t G_f(t-t^\prime) (\boldsymbol{\nabla} {\boldsymbol{v}}_f (t^\prime)+ \boldsymbol{\nabla}{\boldsymbol{v}}_f^{\rm T} (t^\prime)) \mathrm{d} t^\prime, \end{equation}

where $G_f(t)$ is the fluid's shear modulus, and superscript ${\rm T}$ denotes the transpose operation. Similarly, we model the network stress $\boldsymbol {\sigma }_n$ using a general linear VE isotropic CE:

(2.3)\begin{equation} {\boldsymbol{\sigma}}_n = \int_{0}^t [G(t-t^\prime) \left(\boldsymbol{\nabla}{\boldsymbol{v}}_n (t^\prime)+ \boldsymbol{\nabla}{\boldsymbol{v}}_n^{\rm T} (t^\prime) \right) + \lambda(t-t^\prime)\left(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{v}_n(t^\prime)\right)\boldsymbol{I}] \mathrm{d} t^\prime, \end{equation}

where, $G(t)$ and $\lambda (t)$ are the time-dependent first and second Lamé coefficients. For a linear elastic network, these coefficients are related through the Poisson ratio $\nu$: $\lambda ={2G\nu }/{(1-2\nu )}$. Taking the Laplace transform of (2.2) and (2.3) yields $\tilde {\boldsymbol {\sigma }}_f=\tilde {\eta }(s)\,(\boldsymbol {\nabla } \tilde {\boldsymbol {v}}_f+\boldsymbol {\nabla } \tilde {\boldsymbol {v}}^{\rm T}_f)$ and $\tilde {\boldsymbol {\sigma }}_n(s)=\tilde {G}(s)\,(\boldsymbol {\nabla } \tilde {\boldsymbol {v}}_n+\boldsymbol {\nabla } \tilde {\boldsymbol {v}}^{\rm T}_n)+\tilde {\lambda }(s)\,(\boldsymbol {\nabla }\boldsymbol {\cdot } \tilde {\boldsymbol {v}}_n)$, where $\sim$ denotes variables in Laplace space ($s$-space), and $\tilde {\eta }(s)=\mathcal {L}[G_f(t)]$. Substituting these expressions in the Laplace transform of (2.1) gives

(2.4a)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \left( {\tilde{\boldsymbol{v}}}_n \phi + {\tilde{\boldsymbol{v}}}_f (1- \phi) \right) =0 , \end{gather}$$

(2.4b)$$\begin{gather}\tilde{\eta} (s)\,{\nabla}^2 {\tilde{\boldsymbol{v}}}_f + \xi (\tilde{\boldsymbol{v}}_n - \tilde{\boldsymbol{v}}_f ) - (1- \phi )\,\boldsymbol{\nabla} \tilde{p} =\boldsymbol{0} , \end{gather}$$

(2.4c)$$\begin{gather}\tilde{G}(s)\,{\nabla}^2 \tilde{\boldsymbol{v}}_n + (\tilde{\lambda}(s) +\tilde{G}(s)) \boldsymbol{\nabla} (\boldsymbol{\nabla}\boldsymbol{\cdot}\tilde{\boldsymbol{v}}_n ) - \xi (\tilde{\boldsymbol{v}}_n - \tilde{\boldsymbol{v}}_f ) - \phi\, \boldsymbol{\nabla} \tilde{p}=\boldsymbol{0}. \end{gather}$$

Next, we introduce $q:= -\boldsymbol {\nabla }\boldsymbol {\cdot } \tilde {\boldsymbol {v}}_n$. Using (2.4a), we get $\boldsymbol {\nabla } \boldsymbol {\cdot } \tilde {\boldsymbol {v}}_f= -\varepsilon q$, where $\varepsilon = {\phi }/{(\phi -1)}$. We then take the divergence of (2.4b), and (2.4c), and add them to get

(2.5)\begin{equation} {\nabla}^2 \varPhi =0 ,\quad \text{where}\ \varPhi := (\varepsilon \tilde{\eta} + (\tilde{\lambda} +2\tilde{G}) ) q + \tilde{p}. \end{equation}

The general solution of the Laplace equation in spherical coordinates is

(2.6)\begin{equation} \varPhi(r,\theta,\varphi)=\sum_{\ell=0}^\infty \sum_{m={-}\ell}^{m=\ell} \left[D_{\ell ,m}^{{\pm}} \begin{pmatrix} r^{\ell} \\ r^{-\ell -1} \end{pmatrix}\right] Y_{\ell ,m}(\theta,\varphi), \end{equation}

where $Y_{\ell, m}=P_{\ell }^m (\cos \theta )\,{\mathrm {e}}^{{\rm i} m \varphi }$ is the scalar spherical harmonic function of degree $\ell$ and order $m$, $P_{\ell }^m$ is the associated Legendre function of degree $\ell$ and order $m$, $\theta \in [0,{\rm \pi} ]$ is the polar angle, and $\varphi \in [0,2{\rm \pi} )$ is the azimuthal angle; $r^{\ell }$ and $r^{-\ell -1}$ are two linearly independent functions describing the variation of $\varPhi$ in the radial direction. For brevity, we have presented these two functions in the array format $(\ )^{\pm }$, where $D^+$ ($D^-$) are the coefficients associated with the function in the first (second) row, to be determined through imposing BCs. The first row contains the internal solutions (the functions are finite when $r\to 0$ and are unbounded as $r\to \infty$), and the second row contains the external solutions (the functions decay to zero as $r\to \infty$ and are singular as $r\to 0$). This notation is used throughout the paper.

Eliminating the pressure term, after applying the divergence operator to the network and fluid phases, produces a modified Helmholtz equation for $q$:

(2.7)\begin{equation} {\nabla}^2 q - {\alpha}_{\varepsilon}^2 q=0,\quad \text{where}\ {\alpha}_{\varepsilon}^2(s) = \frac{\xi {(1- \varepsilon)}^2}{ \tilde{\lambda}(s) +2\tilde{G}(s) + {\varepsilon}^2 \tilde{\eta} }.\end{equation}

The general solution for $q$, in spherical coordinates, can be expressed as

(2.8)\begin{equation} q(r,\theta,\varphi)=\sum_{\ell=0}^\infty \sum_{m={-}\ell}^{m=\ell}\left[ E_{\ell ,m}^{{\pm}} \begin{pmatrix} {{\mathsf{i}}}_{\ell } ({\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ({\alpha}_{\varepsilon} r) \end{pmatrix} \right] Y_{\ell ,m} (\theta , \varphi), \end{equation}

where ${{\mathsf{i}}}_{\ell } ({\alpha }_{\varepsilon } r)= \sqrt {{{\rm \pi} }/{2 {\alpha }_{\varepsilon } r}}\,I_{\ell +{1}/{2}} ({\alpha }_{\varepsilon } r)$, and ${{\mathsf{k}}}_{\ell } ({\alpha }_{\varepsilon } r)= \sqrt {{2}/{{\rm \pi} {\alpha }_{\varepsilon } r}}\,K_{\ell +{1}/{2}} ({\alpha }_{\varepsilon } r)$, are modified spherical Bessel functions of the first and second kind, respectively (Arfken & Weber Reference Arfken and Weber1999). Taking the general solution of the Laplace equation for $\varPhi$, and the modified Helmholtz equation for $q$, in spherical coordinates, we can express the pressure as

(2.9)\begin{align} \tilde{p} (r,\theta , \varphi)= \sum_{\ell=0}^\infty \sum_{m={-}\ell}^{m=\ell} \left[D_{\ell ,m}^{{\pm}} \begin{pmatrix} r^{\ell} \\ r^{-\ell -1} \end{pmatrix} - ( \varepsilon \tilde{\eta} + \tilde{\lambda} +2\tilde{G} ) \begin{pmatrix} {{\mathsf{i}}}_{\ell } ({\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ({\alpha}_{\varepsilon} r) \end{pmatrix} E_{\ell ,m}^{{\pm}} \right] Y_{\ell ,m} (\theta , \varphi). \end{align}

Next, we introduce ${\boldsymbol {v}}^{+} := \tilde {\eta } \tilde {\boldsymbol {v}}_f + \tilde {G} \tilde {\boldsymbol {v}}_n$ and ${\boldsymbol {v}}^{-} := \tilde {\boldsymbol {v}}_f - \tilde {\boldsymbol {v}}_n$. After summation and subtraction of (2.4b) and (2.4c), and a few lines of algebra, we arrive at the following equations for $\boldsymbol {v}^{+}$ and $\boldsymbol {v}^{-}$:

(2.10a)$$\begin{gather} {\nabla}^2 {\boldsymbol{v}}^{+} - \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{v}}^{+}) = \boldsymbol{\nabla} \varPhi, \end{gather}$$

(2.10b)$$\begin{gather}{\nabla}^2 {\boldsymbol{v}}^{-} + \frac{{\gamma}_{\varepsilon}}{1- \varepsilon}\,\boldsymbol{\nabla} (\boldsymbol{\nabla}\boldsymbol{\cdot} {\boldsymbol{v}}^{-} ) - {\beta}^2 {\boldsymbol{v}}^{-} = \frac{1}{\widetilde{{\eta}}_{\varepsilon}}\,\boldsymbol{\nabla} \varPhi, \end{gather}$$

where

\begin{align*} {\beta }^2 = \xi \left(\frac{1}{\tilde {G}} + \frac{1}{\tilde {\eta }}\right),\quad \frac{1}{{\tilde {\eta }}_{\varepsilon }} = \frac{1}{1-\varepsilon} \left(\frac{1}{\tilde {\eta }} + \frac{\varepsilon }{\tilde {G}}\right) \end{align*}

and

\begin{align*} {\gamma }_{\varepsilon } =\frac{\tilde {\lambda } +\tilde {G}}{\tilde {G}} + \frac{\tilde {\eta } }{1- \varepsilon } \left[ \left( \varepsilon + \frac{\tilde {\lambda } +2\tilde {G}}{\tilde {\eta}} \right) \left(\frac{1}{\tilde {\eta }} + \frac{\varepsilon }{\tilde {G}}\right)\right]. \end{align*}

We note that ${\alpha }_{\varepsilon }^2 = {(1- \varepsilon )}{\beta }^2/{(1- \varepsilon + {\gamma }_{\varepsilon })}$. The general solutions for ${\boldsymbol {v}}^{+}$ and ${\boldsymbol {v}}^{-}$ are:

(2.11)\begin{align} {\boldsymbol{v}}^{+} & = \sum_{\ell=0}^\infty \sum_{m={-}\ell}^{m=\ell} \left\{\begin{pmatrix} r^{\ell} \\ r^{-\ell -1} \end{pmatrix} \left[ A_{\ell ,m}^{{\pm}} \sqrt{\ell (\ell +1)}\,{\boldsymbol{C}}_{\ell ,m} \vphantom{\begin{pmatrix} (\ell +1) {\boldsymbol{P}}_{\ell +1,m} + \sqrt{(\ell +1)(\ell +2)} {\boldsymbol{B}}_{\ell +1,m} \\ - \ell {\boldsymbol{P}}_{\ell -1,m} + \sqrt{\ell (\ell -1)} {\boldsymbol{B}}_{\ell -1,m} \end{pmatrix}}\right.\right. \nonumber\\ &\quad \left.\left.{}+ B_{\ell ,m}^{{\pm}} \begin{pmatrix} (\ell +1) {\boldsymbol{P}}_{\ell +1,m} + \sqrt{(\ell +1)(\ell +2)}\,{\boldsymbol{B}}_{\ell +1,m} \\ - \ell {\boldsymbol{P}}_{\ell -1,m} + \sqrt{\ell (\ell -1)}\,{\boldsymbol{B}}_{\ell -1,m} \end{pmatrix}\right] \right.\nonumber\\ &\quad + \left. D_{\ell ,m}^{{\pm}} \begin{pmatrix} r^{\ell +1} \\ r^{-\ell } \end{pmatrix} \begin{pmatrix} \dfrac{1}{2(2\ell +3)} \\ \dfrac{1}{2(2\ell -1)} \end{pmatrix} \begin{pmatrix} \ell {\boldsymbol{P}}_{\ell ,m} + \dfrac{\ell +3 }{\ell +1} \sqrt{\ell (\ell +1)}\,{\boldsymbol{B}}_{\ell ,m} \\ (\ell +1) {\boldsymbol{P}}_{\ell ,m} + \dfrac{2- \ell}{ \ell} \sqrt{\ell (\ell +1)}\,{\boldsymbol{B}}_{\ell ,m} \end{pmatrix} \right.\nonumber\\ &\quad - \left. \frac{\varepsilon \tilde{\eta} +\tilde{G}}{{\alpha}_{\varepsilon}^2}\, E_{\ell ,m}^{{\pm}} \left[ \frac{\mathrm{d}}{\mathrm{d} r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } ({\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ({\alpha}_{\varepsilon} r) \end{pmatrix} {\boldsymbol{P}}_{\ell ,m}+ \frac{1}{r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } ({\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ({\alpha}_{\varepsilon} r) \end{pmatrix} \sqrt{\ell(\ell +1)}\,{\boldsymbol{B}}_{\ell ,m}\right]\right\}, \end{align}

(2.12)\begin{align} {\boldsymbol{v}}^{-} & = \sum_{\ell=0}^\infty \sum_{m={-}\ell}^{m=\ell} \left\{{A'}_{\ell ,m}^{{\pm}} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix} \sqrt{\ell (\ell +1)}\,{\boldsymbol{C}}_{\ell ,m} \right. \nonumber\\ &\quad \left.{}+ {B '}_{\ell ,m}^{{\pm}} \left[\ell (\ell +1) {\boldsymbol{P}}_{\ell ,m}\,\frac{1}{\beta r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix} + \sqrt{\ell (\ell +1)}\,{\boldsymbol{B}}_{\ell ,m} \right.\right. \nonumber\\ &\quad \left.{}\times \left(\frac{\mathrm{d}}{\mathrm{d} (\beta r)} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix} + \frac{1}{\beta r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix}\right)\right] \nonumber\\ &\quad - \frac{1}{ {\tilde{\eta}}_{\varepsilon} {\beta}^2}\,{D}_{\ell ,m}^{{\pm}} \begin{pmatrix} r^{\ell -1} \\ r^{-\ell -2} \end{pmatrix} \begin{pmatrix} \ell {\boldsymbol{P}}_{\ell ,m} + \sqrt{\ell (\ell +1)}\,{\boldsymbol{B}}_{\ell ,m} \\ -(\ell +1) {\boldsymbol{P}}_{\ell ,m} + \sqrt{\ell (\ell +1)}\,{\boldsymbol{B}}_{\ell ,m} \end{pmatrix} \nonumber\\ &\quad \left.{}-\frac{ \varepsilon -1 }{ {\alpha}_{\varepsilon}^2 }\,{E}_{\ell ,m}^{{\pm}} \left[\frac{\mathrm{d}}{\mathrm{d} r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } ( {\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ( {\alpha}_{\varepsilon} r) \end{pmatrix} {\boldsymbol{P}}_{\ell ,m} + \frac{1}{r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } ( {\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ( {\alpha}_{\varepsilon} r) \end{pmatrix} \sqrt{\ell (\ell +1)}\,{\boldsymbol{B}}_{\ell ,m}\right]\right\}, \end{align}

where

(2.13a)$$\begin{gather} {\boldsymbol{P}}_{\ell ,m} = Y_{\ell ,m} (\theta ,\varphi)\,\hat{\boldsymbol{r}} , \end{gather}$$

(2.13b)$$\begin{gather}{\boldsymbol{B}}_{\ell ,m} =\frac{1}{\sqrt{\ell (\ell +1)}} \left[ \frac{\partial }{\partial \theta} \hat{\boldsymbol{\theta}} + \frac{1}{\sin \theta}\,\frac{\partial }{\partial \varphi} \widehat{\boldsymbol{\varphi}} \right] Y_{\ell,m} (\theta , \varphi) , \end{gather}$$

(2.13c)$$\begin{gather}{\boldsymbol{C}}_{\ell ,m} = \frac{1}{\sqrt{\ell (\ell +1)}} \left[\frac{1}{\sin \theta}\,\frac{\partial }{\partial \varphi} \hat{\boldsymbol{\theta}} - \frac{\partial }{\partial \theta} \widehat{\boldsymbol{\varphi}} \right] Y_{\ell,m} (\theta , \varphi) \end{gather}$$

are scaloidal–poloidal–toroidal representations of mutually orthogonal vector spherical harmonics that are used commonly to express the solutions to vector Laplace and Helmholtz equations in spherical coordinates (Morse & Feshbach Reference Morse and Feshbach1953), and $A_{\ell,m}^{\pm }$, $B_{\ell,m}^{\pm }$, ${A'}_{\ell,m}^{\pm }$, ${B'}_{\ell,m}^{\pm }$, $D_{\ell,m}^{\pm }$ and $E_{\ell,m}^{\pm }$ are 12 constants to be determined from BCs; the solution details are given in Appendix A. Substituting these expressions in relations $\tilde {\boldsymbol {v}}_f = ({\boldsymbol {v}}^{+} + \tilde {G} {\boldsymbol {v}}^{-} )/(\tilde {\eta } + \tilde {G})$ and $\tilde {\boldsymbol {v}}_n = ({\boldsymbol {v}}^{+} - \eta {\boldsymbol {v}}^{-} )/(\tilde {\eta } + \tilde {G})$ gives the fluid and network velocity (or displacement) fields in $s$-space, which we do not produce here for brevity.

2.1. Axisymmetric solutions

The general solutions for $\tilde {\boldsymbol {v}}_f$ and $\tilde {\boldsymbol {v}}_n$ can be used to obtain analytical expressions for a variety of problems involving spherical geometries. Several applications of Lamb's general solution for Stokes flow may be found in Happel & Brenner (Reference Happel and Brenner2012) and Kim & Karrila (Reference Kim and Karrila2013); the analytical solutions of poroelastic equations using Biot's general solution, which neglects the fluid shear forces, are surveyed in Detournay & Cheng (Reference Detournay and Cheng1993) and Cheng (Reference Cheng2016). The majority of these solutions involve axisymmetric geometries, i.e. $\partial _{\varphi } v_r = \partial _{\varphi } v_{\theta } =0$, and ${v}_\varphi =0$. This is achieved when the summation index is $m=0$ in (2.11) and (2.12), and ${A}_{\ell m}^{\pm }={A'}_{\ell m}^{\pm }=0$. The axisymmetric forms of the fluid and network velocities in Laplace space simplify to

(2.14a)\begin{align} \tilde{v}_{r,f} &= \frac{1}{\tilde{\eta} +\tilde{G}} \left[\sum_{\ell =0}^{\infty} \left\{\begin{pmatrix} {B}_{\ell -1}^{+}\\{B}_{\ell +1}^{-} \end{pmatrix} \begin{pmatrix} \ell r^{\ell -1} \\ -(\ell +1) r^{-\ell -2} \end{pmatrix} + {B'}_{\ell}^{{\pm}}\, \frac{\ell (\ell +1) \tilde{G}}{\beta r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix} \vphantom{\begin{pmatrix} \dfrac{\ell +3}{2(\ell +1)(2\ell +3)} r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \\ \dfrac{2- \ell}{2 \ell (2\ell -1)} r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \end{pmatrix}}\right.\right.\nonumber\\ &+ D_{\ell }^{{\pm}} \begin{pmatrix} r^{\ell -1} \\ r^{-\ell -2} \end{pmatrix} \begin{pmatrix} \dfrac{\ell}{2(2\ell +3)}\,r^2 - \dfrac{G}{{\tilde{\eta}}_{\varepsilon} {\beta}^2}\,\ell \\ \dfrac{\ell +1}{2(2\ell -1)}\,r^2 + \dfrac{\tilde{G}}{ {\tilde{\eta}}_{\varepsilon} {\beta}^2}\,(\ell +1) \end{pmatrix} \left.\left.\vphantom{\begin{pmatrix} \dfrac{\ell +3}{2(\ell +1)(2\ell +3)} r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \\ \dfrac{2- \ell}{2 \ell (2\ell -1)} r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \end{pmatrix}} - \frac{\varepsilon (\tilde{\eta} +\tilde{G})}{{\alpha}_{\varepsilon}^2}\,{E}_{\ell}^{{\pm}}\, \frac{\mathrm{d}}{\mathrm{d}r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } ( {\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ( {\alpha}_{\varepsilon} r) \end{pmatrix} \right\} P_{\ell} (\cos\theta) \right], \end{align}

(2.14b)\begin{align} \tilde{v}_{\theta ,f} &= \frac{1}{\tilde{\eta} +\tilde{G}} \left[\sum_{\ell =1 }^{\infty} \left\{ \begin{pmatrix} {B}_{\ell -1}^{+} \\ {B}_{\ell +1}^{-} \end{pmatrix} \begin{pmatrix} r^{\ell -1} \\ r^{-\ell -2} \end{pmatrix} + {B'}_{\ell}^{{\pm}}\,\frac{\tilde{G}}{\beta } \left[\frac{\mathrm{d}}{\mathrm{d}r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix} + \frac{1}{r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix} \right] \vphantom{\begin{pmatrix} \dfrac{\ell +3}{2(\ell +1)(2\ell +3)} r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \\ \dfrac{2- \ell}{2 \ell (2\ell -1)} r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \end{pmatrix}}\right.\right.\nonumber\\ &\quad +\left.\left. D_{\ell }^{{\pm}} \begin{pmatrix} r^{\ell -1} \\ r^{-\ell -2} \end{pmatrix} \begin{pmatrix} \dfrac{\ell +3}{2(\ell +1)(2\ell +3)}\,r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \\ \dfrac{2- \ell}{2 \ell (2\ell -1)}\,r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \end{pmatrix} - \frac{\varepsilon (\tilde{\eta} +\tilde{G})}{{\alpha}_{\varepsilon}^2}\,{E}_{\ell}^{{\pm}}\,\frac{1}{r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } ({\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ( {\alpha}_{\varepsilon} r) \end{pmatrix}\right\} \frac{\mathrm{d}}{\mathrm{d}\theta} P_{\ell} (\cos\theta)\right], \end{align}

(2.15a)\begin{align} \tilde{v}_{r,n} &= \frac{1}{\tilde{\eta} + \tilde{G}} \left[\sum_{\ell =0}^{\infty} \left\{ \begin{pmatrix} {B}_{\ell -1}^{+}\\ {B}_{\ell +1}^{-} \end{pmatrix} \begin{pmatrix} \ell r^{\ell -1} \\ -(\ell +1) r^{-\ell -2} \end{pmatrix} - {B'}_{\ell}^{{\pm}}\,\frac{\ell (\ell +1) \tilde{\eta}}{\beta r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix}\right.\right.\nonumber\\ &\quad + \left.\left. D_{\ell }^{{\pm}} \begin{pmatrix} r^{\ell -1} \\ r^{-\ell -2} \end{pmatrix} \begin{pmatrix} \dfrac{\ell}{2(2\ell +3)}\,r^2 + \dfrac{\tilde{\eta}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2}\,\ell \\ \dfrac{\ell +1}{2(2\ell -1)}\,r^2 - \dfrac{\tilde{\eta}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2}\,(\ell +1) \end{pmatrix} - \frac{ \tilde{\eta} + \tilde{G} }{{\alpha}_{\varepsilon}^2}\,{E}_{\ell}^{{\pm}}\, \frac{\mathrm{d}}{\mathrm{d}r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } ( {\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ( {\alpha}_{\varepsilon} r) \end{pmatrix} \right\} P_{\ell} (\cos\theta) \right], \end{align}

(2.15b)\begin{align} \tilde{v}_{\theta ,n} &= \frac{1}{\tilde{\eta} +\tilde{G}} \left[ \sum_{\ell =1 }^{\infty} \left\{ \begin{pmatrix} {B}_{\ell -1}^{+} \\ {B}_{\ell +1}^{-} \end{pmatrix} \begin{pmatrix} r^{\ell -1} \\ r^{-\ell -2} \end{pmatrix} - {B'}_{\ell}^{{\pm}}\,\frac{\tilde{\eta}}{\beta } \left[\frac{\mathrm{d}}{\mathrm{d}r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix} + \frac{1}{r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } (\beta r) \\ {{\mathsf{k}}}_{\ell } (\beta r) \end{pmatrix} \right] \vphantom{\begin{pmatrix} \dfrac{\ell +3}{2(\ell +1)(2\ell +3)} r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \\ \dfrac{2- \ell}{2 \ell (2\ell -1)} r^2 - \dfrac{\tilde{G}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \end{pmatrix}}\right.\right.\nonumber\\ &\quad +\left.\left. D_{\ell }^{{\pm}} \begin{pmatrix} r^{\ell -1} \\ r^{-\ell -2} \end{pmatrix} \begin{pmatrix} \dfrac{\ell +3}{2(\ell +1)(2\ell +3)}\,r^2 + \dfrac{\tilde{\eta}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \\ \dfrac{2- \ell}{2 \ell (2\ell -1)}\,r^2 + \dfrac{\tilde{\eta}}{{\tilde{\eta}}_{\varepsilon} {\beta}^2} \end{pmatrix} - \frac{ \tilde{\eta} + \tilde{G} }{{\alpha}_{\varepsilon}^2}{E}_{\ell}^{{\pm}}\,\frac{1}{r} \begin{pmatrix} {{\mathsf{i}}}_{\ell } ( {\alpha}_{\varepsilon} r) \\ {{\mathsf{k}}}_{\ell } ( {\alpha}_{\varepsilon} r) \end{pmatrix} \right\} \frac{\mathrm{d}}{\mathrm{d}\theta} P_{\ell} (\cos\theta)\right]. \end{align}

We note that fluid velocity dependency on the network compressibility (through ${\alpha }_{\varepsilon }$) is proportional to $\phi \sim \varepsilon \ll 1$ and thus negligible in most physiologically relevant conditions. In the next section, we use the axisymmetric solutions to study the dynamics of a spherical bead moving within a PVE medium, and provide an exact solution of the equations subject to appropriate BCs. Moreover, we discuss the differences between the dynamics of the sphere moving in VE and PVE materials in the context of single-particle and two-point MR.

3. A rigid sphere moving a network that is permeated by a fluid

Consider a rigid sphere moving with time-dependent velocity $U (t)\, \hat {\boldsymbol {z}}$ in an unbounded PVE medium. We assume that the size of the sphere is larger than the typical mesh size of the background network, and assume no-slip BCs for both the network and fluid velocity fields:

(3.1a,b)\begin{equation} \text{at}\ r=a:\quad \tilde{v}_{r,f}=\tilde{v}_{r,n}=\tilde{U}(s)\cos \theta,\quad \tilde{v}_{\theta,f}=\tilde{v}_{\theta,n} ={-}\tilde{U}(s)\sin \theta. \end{equation}

Note that these BCs are accurate only when the interface is between an elastic (or rigid) domain and a PVE domain. The BCs become more involved at fluid–PVE and PVE–PVE interfaces (Feng & Young Reference Feng and Young2020). Since the velocity field must decay to zero at infinity, and functions corresponding to internal solutions are all unbounded as $r\to \infty$, their corresponding coefficients ($+$ superscript) are identically zero.

Applying no-slip BCs (3.1a,b) to the axisymmetric solutions for the fluid and network, we get

(3.2a)\begin{align} \tilde{v}_{r,f} &= \frac{a}{{\beta}^2 r^3}\,\frac{\tilde{U}(s)}{{\varDelta}_{\varepsilon}} \left[ g_1 + (1-\varepsilon) {\tau}_{\varepsilon} (3 r^2 - a^2) {\beta}^2 g_2 -3\, \frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\,(1-\varepsilon) \,{\rm e}^{-\beta (r-a)} (1+ \beta r) \right. \nonumber\\ &\quad \left.\vphantom{\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}} -\varepsilon g_3 + 3\varepsilon (1+\tau) \exp({- {\alpha}_{\varepsilon} (r-a)}) (2+ 2{\alpha}_{\varepsilon} r + {\alpha}_{\varepsilon}^2 r^2) \right] \cos\theta, \end{align}

(3.2b)\begin{align} \tilde{v}_{\theta ,f} &= \frac{a}{2{\beta}^2 r^3}\,\frac{\tilde{U}(s)}{{\varDelta}_{\varepsilon}} \left[g_1 - (1-\varepsilon) {\tau}_{\varepsilon} (3 r^2 + a^2) {\beta}^2 g_2 \vphantom{\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}}\right. \nonumber\\ &\quad -3\,\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\,(1-\varepsilon) \,{\rm e}^{-\beta (r-a)} (1+ \beta r + {\beta}^2 r^2) \nonumber\\ &\quad\left.\vphantom{\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}} -\varepsilon g_3 + 6\varepsilon (1+\tau) \exp({- {\alpha}_{\varepsilon} (r-a)}) (1+ {\alpha}_{\varepsilon} r ) \right] \sin\theta, \\ \tilde{v}_{r,n} &= \frac{a}{{\beta}^2 r^3}\,\frac{\tilde{U}(s)}{{\varDelta}_{\varepsilon}} \left[{-}g_4 + (1-\varepsilon) {\tau}_{\varepsilon} (3 r^2 - a^2) {\beta}^2 g_2 +3\, \frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\,\tau (1-\varepsilon) \,{\rm e}^{-\beta (r-a)} (1+ \beta r) \right. \nonumber\\ &\quad\left.\vphantom{\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}} +\varepsilon \tau g_1 + 3 (1+\tau) \exp({- {\alpha}_{\varepsilon} (r-a)}) (2+ 2{\alpha}_{\varepsilon} r + {\alpha}_{\varepsilon}^2 r^2) \right] \cos\theta, \\ \tilde{v}_{\theta ,n} &= \frac{a}{2{\beta}^2 r^3}\,\frac{\tilde{U}(s)}{{\varDelta}_{\varepsilon}} \left[{-}g_4 - (1-\varepsilon) {\tau}_{\varepsilon} (3 r^2 + a^2) {\beta}^2 g_2 \vphantom{\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}}\right. \nonumber\\ &\quad +3\,\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\,\tau (1-\varepsilon) \,{\rm e}^{-\beta (r-a)} (1+ \beta r + {\beta}^2 r^2) \nonumber\\ &\quad \left.\vphantom{\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}} +\varepsilon \tau g_1 + 6 (1+\tau) \exp({- {\alpha}_{\varepsilon} (r-a)}) (1+ {\alpha}_{\varepsilon} r ) \right] \sin\theta, \end{align}

where $\tau (s) = {\tilde {\eta }(s)}/{\tilde {G} (s)}$ is the shear relaxation time, $\tilde {\eta }_{\varepsilon } = {(1-\varepsilon ) \tilde {\eta }}/{(1+ \varepsilon \,\tau (s))}$, and ${\tau }_{\varepsilon } (s)= {\tilde {\eta }_{\varepsilon }}/{\tilde {G}(s)}$. Defining the inverse of permeability as ${\beta }_{\circ }^2 = {\xi }/{\tilde {\eta }}$, we have

\[ {\beta }^2 = {\beta }_{\circ }^{2} (1+\tau (s) ), {\alpha }_{\varepsilon }^2 = {\beta }_{\circ }^{2} \frac{\tau (s)(1-\varepsilon )^2}{\frac{2(1-\nu )}{1-2\nu} + {\varepsilon }^2\,\tau (s)}. \]

The remaining coefficients that appear in (3.2) are given by the expressions

(3.3a)$$\begin{gather} {\varDelta}_{\varepsilon} = 2 {\tau}_{\varepsilon} (1-\varepsilon) (1+ a {\alpha}_{\varepsilon}) + \frac{{\alpha}_{\varepsilon}^2}{{\beta}^2} (1+ \varepsilon \tau + {\tau}_{\varepsilon} (1-\varepsilon) (1+ a\beta + a^2 {\beta}^2)), \end{gather}$$

(3.3b)$$\begin{gather}g_1 = \frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\,(3+3a \beta + a^2 {\beta}^2),\quad g_2 = 1+ a {\alpha}_{\varepsilon} + \frac{1}{2}\,\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\, (1+ a\beta + a^2 {\beta}^2), \end{gather}$$

(3.3c)$$\begin{gather}g_3 = 6(1+\tau) (1+ a {\alpha}_{\varepsilon}) + 2 \tau a^2 {\alpha}_{\varepsilon}^2 + 3\, \frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\,(1+ a\beta + a^2 {\beta}^2), \end{gather}$$

(3.3d)$$\begin{gather}g_4 = 6(1+\tau) (1+ a {\alpha}_{\varepsilon}) + 2 a^2{\alpha}_{\varepsilon}^2 + 3 \tau\,\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\,(1+ a\beta + a^2 {\beta}^2). \end{gather}$$

Substituting for coefficients in (2.9), the pressure is given by

(3.4)\begin{align} \tilde{p} &= \frac{3 a (\tilde{\eta} +\tilde{G})}{2 r^2}\, \frac{\tilde{U}(s)}{{\varDelta}_{\varepsilon}} \left[2 {\tau}_{\varepsilon} (1-\varepsilon) (1+ a {\alpha}_{\varepsilon}) + {\tau}_{\varepsilon}\,\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\, (1- \varepsilon) (1+a\beta + a^2 {\beta}^2) \right. \nonumber\\ &\quad \left. {}-2\,\frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\left(\varepsilon \tau + \frac{\tilde{\lambda}}{\tilde{G}} +2 \right) \exp({-{\alpha}_{\varepsilon} (r-a)}) (1+ {\alpha}_{\varepsilon} r) \right] \cos\theta. \end{align}

The fluid and network stress components, including the traceless and isotropic parts (see Appendix B for explicit expressions of the network and fluid stresses), are

(3.5a–d)\begin{gather} \tilde{\sigma}_{rr}^{f} ={-}\tilde{p} (1-\phi)+ 2 \tilde{\eta}\, \frac{\partial \tilde{v}_{r,f}}{\partial r},\\ \tilde{\sigma}_{rr}^{n} = (\tilde{\lambda} + 2\tilde{G})\,\frac{\partial \tilde{v}_{r,n} }{\partial r} + \frac{\tilde{\lambda}}{r} \left( 2 \tilde{v}_{r,n} + \cot\theta\,\tilde{v}_{ \theta ,n} + \frac{\partial}{\partial \theta} \tilde{v}_{ \theta ,n} \right) -\phi \tilde{p}, \\ \tilde{\sigma}_{r\theta}^{f} = \tilde{\eta} \left( \frac{1}{r}\, \frac{\partial \tilde{v}_{r,f}}{\partial \theta} + \frac{\partial \tilde{v}_{\theta , f}}{\partial r} - \frac{\tilde{v}_{\theta , f}}{r}\right),\\ \tilde{\sigma}_{r\theta}^{n}= \tilde{G} \left( \frac{1}{r}\,\frac{\partial \tilde{v}_{r,n}}{\partial\theta} + \frac{\partial \tilde{v}_{ \theta ,n}}{\partial r} - \frac{\tilde{v}_{ \theta ,n}}{r} \right). \end{gather}

Evaluating these stresses at $r=a$ and integrating over the sphere surface gives the following relations for the total force from the fluid and network phases:

(3.6)\begin{align} \tilde{\boldsymbol{F}}(s) &= \int_{0}^{2{\rm \pi}} \int_{0}^{\rm \pi}((\tilde{\sigma}_{rr}^{f} + \tilde{\sigma}_{rr}^{n}) \hat{\boldsymbol{r}} + ( \tilde{\sigma}_{r\theta}^{f} + \tilde{\sigma}_{r\theta}^{n} ) \hat{\boldsymbol{\theta}}) a^2 \sin\theta \,\mathrm{d}\theta\,\mathrm{d}\varphi =\tilde{R}(s)\,\tilde{U}(s)\,\hat{\boldsymbol{z}}, \end{align}

(3.7)\begin{align} \tilde{R}(s) &= 6{\rm \pi}\,\tilde{G}(s)\,a \frac{\left(1+\tau\right)\tau}{{\varDelta}_{\varepsilon}} \left[\frac{\tilde{\eta}_{\varepsilon}}{\tilde\eta}\,(1-\varepsilon) \left(2(1+a {\alpha}_{\varepsilon} ) + \frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\,(1+ a\beta + a^2 {\beta}^2)\right) \right.\nonumber\\ &\quad \left.{}+ \frac{2}{3}\,\varepsilon\, \frac{{\alpha}_{\varepsilon}^2}{{\beta}^2}\,(1+ a {\alpha}_{\varepsilon})\right], \end{align}

where $\tilde {R}(s)$ is the response function that we set out to find. Using the fluctuation dissipation theorem (FDT), the measurable mean squared displacement (MSD) can be related to the response function using $\langle {\rm \Delta} \tilde {\boldsymbol {r}}^2 \rangle =6k_{B}T\,\tilde {M}(s)/s^2$ (Squires & Mason Reference Squires and Mason2010), where $\langle \,\rangle$ denotes ensemble average, ${\rm \Delta} \tilde {\boldsymbol {r}}^2$ is the MSD of the spherical probe, $\tilde {M}(s)=\tilde {R}^{-1}(s)$ is the sphere's isotropic mobility in Laplace space, and $k_{B} T$ is Boltzmann thermal energy. This relation is used to compute $\tilde {R}(s)$, which in theory can be used to compute $\tilde {G}(s)$, $\tilde {\eta }(s)$, $\nu$ and $\beta _\circ$ in single-particle tracking MR.

In two-point MR, the positional cross-correlations of two fluctuating probes are used to measure the properties of the medium separating them. We can relate these movements to the medium's mechanical properties by applying the FDT, which yields $\mathcal {L}\langle {\rm \Delta} \boldsymbol {x}_1(0)\,{\rm \Delta} \boldsymbol {x}_2(t)\rangle =6 k_{B} T\,\tilde {\boldsymbol {M}}_{21}(s)/s^2$ (Squires & Mason Reference Squires and Mason2010), where $\mathcal {L}$ is the Laplace transform operator, and $\tilde {\boldsymbol {M}}_{21}(s)$ is the pair mobility tensor that computes the velocity of particle $2$ due to a force on particle $1$: $\tilde {\boldsymbol {U}}_2(\boldsymbol {r},s)=\tilde {\boldsymbol {M}}_{21} (\boldsymbol {r},s)\boldsymbol {\cdot }\tilde {\boldsymbol {F}}_1(s)$, where $\boldsymbol {r}$ is the pair separation vector. For VE materials, and when the particles are separated by a large distance ($r/a\gg 1$), the motion of particle 2 becomes independent of its size and asymptotes to the velocity of the surrounding VE fluid: $\tilde {\boldsymbol {U}}_2(\boldsymbol {r},s)\approx \tilde {\boldsymbol {v}}(\boldsymbol {r},s)$.

We use a similar approximation here and assume that particle 2 moves with the network velocity at its centre: $\boldsymbol {U}_2(\boldsymbol {r})\approx \boldsymbol {v}_n(\boldsymbol {r})$ when $r/a\gg 1$. The pair mobility tensor can be expressed in the general form $\tilde {\boldsymbol {M}}_{21}=\widetilde {{M}}_{21}^\parallel \hat {\boldsymbol {r}}\hat {\boldsymbol {r}}+\widetilde {{M}}_{21}^\perp (\boldsymbol {I}-\hat {\boldsymbol {r}}\hat {\boldsymbol {r}})$, where $\parallel$ and $\perp$ denote parallel and perpendicular to the direction of applied force, respectively. Using $\tilde {U}(s)=\tilde {M}(s)\,\tilde {F}(s)$ in (3.2), and the above expression for the mobility tensor, we arrive at the expressions

(3.8a,b)\begin{equation} \tilde{M}_{21}^\parallel{=}\tilde{M}(s)\,\frac{a}{\beta^2 r^3 \varDelta_\epsilon}\,[\cdots]_{r,n},\quad \tilde{M}_{21}^\perp{=}\tilde{M}(s)\,\frac{a}{2 \beta^2 r^3 \varDelta_\epsilon}\,[\cdots]_{\theta,n}, \end{equation}

where $[\cdots ]_{r,n}$ and $[\cdots ]_{\theta,n}$ are the terms between brackets in (3.2c) and (3.2d). These expressions provide the mathematical framework for analysing two-point MR results.

3.1. Results

In this section, we consider the example of a spherical probe moving under constant force ${\boldsymbol {F}}_{\circ }$, as an analogue to the active MR experiment, within a PVE material composed of a linear elastic network, with shear modulus $G$ and Poisson ratio $\nu$, and a Newtonian fluid of viscosity $\eta$, making $\tau (s)=s\eta /G$; see figure 1(a). The probe's velocity and the induced displacements in the network and fluid phases can be computed in Laplace space using (3.6) and (3.2), respectively. We then use Fourier–Euler summation (Abate & Whitt Reference Abate and Whitt1992) to invert numerically the results from $s$-space to time-space. The limiting values of all these quantities at $t=0$ and $t\to \infty$ can be computed analytically using the following limits in $s$-space: $f(t=0^+) =\lim _{s\to \infty } s\,\tilde {f}(s)$ and $f(t\to \infty ) =\lim _{s\to 0} s\,\tilde {f}(s)$ (see Appendix C). Figure 1(b) shows the variations of $1-x/x_\infty$ versus $t/\tau$ for different values of $\beta _0>1$, when $\phi \to 0$ and $\nu =0.3$ (length is non-dimensionalized with probe radius). Here, $x$ is the probe's position, ${x}_{\infty }=({{F_{\circ }}}/{6{\rm \pi} G a}) ({(5-6\nu )}/{4(1-\nu )})$ is the long-time position of the probe, and $\tau =\eta /G$ is the time scale of shear deformations. Note that the curves in the main figures produce identical results for different values of $\tau$ when time is made dimensionless as $t/\tau$. The predictions of the VE model ($1-x_{{VE}}/{x_\infty }=\textrm {e}^{-t/\tau }$) are displayed as a dashed line. As can be seen, the relaxation at early times is independent of permeability $\beta _{\circ }$, and is controlled by $\tau$. This is followed by a slow relaxation dynamics, with the relaxation time appearing to increase with permeability $\beta _\circ$. Note that $1-x/x_\infty \to 0$ as $t \to \infty$.

Figure 1. (a) A schematic representation of a rigid spherical probe moving under a constant force in a PVE medium. The sphere moves from ${{x}}_{\circ }$ at early times to ${{x}}_{\infty }$ at long times. Plots of the sphere relaxation $1-x(t)/x_\infty$ versus time are shown for: (b) different values of permeability, when $\phi \to 0$ and $\nu =0.3$; (c) different values of the Poisson ratio, when $\phi \to 0$ and ${\beta }_{\circ } =2$; and (d) different values of the network volume fraction, when $\nu =0.3$ and ${\beta }_{\circ } =2$. Each inset figure represents the results of its associated main figure, as a function of rescaled time $t/\tau _D$, where $\tau _D=a^2/D_q$ is the diffusion time scale of the compressibility field for distance $a$, where $D_q=\tau ^{-1}\alpha _\epsilon ^{-2}$ is the diffusion coefficient of the network compressibility.

Next, to study the effect of network compressibility, we compute the relaxation dynamics for different values of the Poisson ratio $\nu$ and a constant permeability $\beta _\circ =2$; see figure 1(c). As can be seen, the extent of deformation that follows a slow relaxation dynamics is reduced with increasing $\nu$. For $\nu =0.49$, which correspond to a nearly incompressible network, the relaxation is given almost entirely by the VE model. This observation suggests that the slow relaxation is induced by network compressibility. An inspection of (2.7) shows that for a linear elastic network and a Newtonian fluid, it is a diffusion equation, with the diffusion coefficient given by $D_q=\tau ^{-1} \alpha _\epsilon ^{-2}$. The same equation and concept appear in Biot poroelasticity (Detournay & Cheng Reference Detournay and Cheng1993; Doi Reference Doi2009; Cheng Reference Cheng2016), where $D_q$ is sometimes referred to as generalized consolidation coefficient in the soil mechanics literature.

Following this observation, we rescale the time axis of figure 1(b) with $t/\tau _D$, where $\tau _D=a^2/D_q = a^2 {\alpha }_{\varepsilon }^2 \tau$ is the characteristic time scale for the diffusion of network compressibility. Upon this rescaling, all the results collapse to a single line at longer time scales (see the inset of figure 1b), which is consistent with the idea that the slow relaxation is determined by the diffusion of network compressibility and the fluid pressure induced by that. This rescaling did not result in the collapse of the data for different values of $\nu$ (see the inset of figure 1c), suggesting a more complex dependency of the dynamics on $\nu$. To further clarify the role of shear viscosity, in Appendix D, we include the bead's relaxation dynamics, where we vary only the shear viscosity and keep the rest of the parameters unchanged. As we explained earlier, decreasing the viscosity leads to faster initial relaxation and does not affect the long-time relaxation as long as $\xi$ remains unchanged.

In figure 1(d), we study the effect of volume fraction of the network phase on probe dynamics in the range $0 < \phi \le 0.3$ for the choice $\beta _\circ =2$ and $\nu =0.3$. The value of $\phi$ appears not to change the qualitative behaviour and has only a minor quantitative effect on the overall dynamics. We observe the same weak dependencies for the displacement fields. As such, in the remainder of this paper, we present only the results for $\phi \to 0$, which results in $\varepsilon =0$, $\eta _\varepsilon =\eta$, $\tau _\varepsilon =\tau$ and $\alpha _\varepsilon =\alpha$.

Note that even though $\phi \to 0$, the fluid permeability can still be largely affected by the cytoskeletal network. To see this, consider the expression for permeability of a dilute fibrous network in a Newtonian fluid: $(\beta _\circ a_f)^{2}={8\phi }/{(-\ln \phi -1.48)}$ (Sangani & Acrivos Reference Sangani and Acrivos1982), where $a_f$ is the radius of the fibre. Let us consider microtubules that have the largest radius among the cytoskeletal filaments: $a_f \approx 12 \,\textrm {nm}$. Rescaling this expression to the probe's radius – which is the scale of interest here – gives $({\beta _\circ } a)^2= ({a}/{a_f})^2({8\phi }/{(-\ln \phi -1.48)})$. Assuming a spherical bead of radius $a\sim 1 \,\mathrm {\mu }\textrm {m}$, we see that even for very small volume fractions, we have the likely condition $\phi (a/a_f)^2 \gg 1$, which results in $\beta _\circ a \gg 1$ and a substantial reduction in fluid permeability by the network in the probe's length scale.

Based on these findings, we provide a recipe for determining the constitutive parameters in an active MR set-up such as the one used here, when the fluid is Newtonian. If the fluid is also VE, then the probe motion alone can determine only the ratio $\tilde {\eta }(s)/\tilde {G}(s)$ and not the individual terms, and displacement fields will be needed to disentangle the time scales associated with each of the fluid and the network. At very early times, the network is hardly deformed and the probe velocity is determined by the fluid drag force: $U (t\to 0)={F}_\circ /(6{\rm \pi} \eta a)$. Measuring $U (t=0)$ can thus be used to compute $\eta$. At later intermediate times the dynamics is controlled by relaxation of shear modes, allowing us to determine $\eta /G(t)$ and thus $G(t)$. Having $\eta$ and $G(t)$, one can compute $\nu$ by using the steady-state displacement of the probe:

\[ {x}_{\infty }=\frac{{F}_{\circ }}{6{\rm \pi} G(t\to \infty ) a} \frac{5-6\nu }{4(1-\nu )}. \]

Finally, $\beta _{\circ }$ can be determined by fitting (3.6) to the displacement versus time at longer time scales.

We now focus on the network displacement versus time at different, but all large, distances from the probe. Figure 2(a) shows displacement relaxation along the direction of the applied force, $1-u_{n}^{\parallel }(r)/u_{n}^{\parallel,\infty }$, at different distances from the probe, for $\beta _\circ =2$, $\nu =0.3$ and $\phi \to 0$ (as for the rest of the results). Similar to the probe's dynamics, the relaxation at early times is determined by the shear modes given by the VE model (dashed line). For large $r/a$, the time-dependent $1-u_{n}^{\parallel }/u_{n}^{\parallel,\infty }$ exhibits a slow relaxation dynamics that varies with separation distance and exhibits a local minimum (a minimum for $u_{n}^{\parallel }$) followed by what appears to be a local maximum, noting that $1-u_{n}^{\parallel }/u_{n}^{\parallel,\infty } \to 0$ as $t\to \infty$.

Figure 2. The relaxation of network displacement fields in response to the moving probe in directions parallel (a) and perpendicular (b) to the applied force, and at different distances from the sphere, $1-u^{(\parallel,\perp )}_n/u^{(\parallel,\perp ),\infty }_n$. The inset figures show the same results as the associated figures when the time axis is scaled as $t/{\tau }_{D}^r$, where $\tau _D^r=r^2 {\alpha }^2 \tau$ is the time scale for the network compressibility field to diffuse distance $r$ away from the probe. Here, we have taken ${\beta }_{\circ }=2$, $\phi \to 0$ and $\nu =0.3$.

Figure 2(b) shows the displacement relaxation in the direction perpendicular to the applied force, $u_{n}^{\perp }$, versus time at different values of $r$. Unlike the displacements in the parallel direction, we do not observe any local optimum in time. Aside from this difference, the relaxation dynamics follows closely the behaviour that we observed in the parallel direction.

The predicted change in relaxation dynamics with $r$ is in line with previous experimental results (Rosenbluth et al. Reference Rosenbluth, Crow, Shaevitz and Fletcher2008), where atomic force microscopy (AFM) was used to deform the cell cortex, and the displacements of probe particles were measured at different distances from the AFM tip. Similar to our predictions, the experiments show a fast initial relaxation that is independent of distance from the AFM tip, followed by a significantly slower relaxation dynamics that slows down further with increasing distance from the tip.

To test if the slow relaxation is again controlled by the diffusion of the stress associated with network compressibility, we scale the time axis with the diffusion time scale for travelling the distance $r$: $\tau _D^{r}=r^2/D_q = r^2 {\alpha }^2 \tau$. Again, we observe a nice collapse of the plots at longer times, which further confirms that slow relaxation dynamics is determined by the diffusion of network compressibility and its associated stresses.

In Appendix D, we present the results of pressure relaxation dynamics induced by the probe's motion under a constant force. In figure 7, we show the relaxation dynamics of the pressure on the surface of the probe and parallel to the applied force ($p(r=a,\theta =0)$), for different values of permeability ${\beta }_{\circ }$ and Poisson ratio $\nu$ of the PVE material. As can be seen, in all cases the pressure relaxation is controlled entirely by $\tau _D$ and is independent of shear relaxation $\tau$. We observe the same behaviour when studying the pressure relaxation at different probe distances (different values of $r/a$); see figure 7. This behaviour is in agreement with the previous two-dimensional immersed boundary simulations (Strychalski & Guy Reference Strychalski and Guy2016). (We note that a recent three-dimensional immersed boundary simulation of cell blebbing by Strychalski (Reference Strychalski2021) observes a sub-diffusive pressure relaxation.)

In figures 4(a,b) of Appendix D, we present the relaxation dynamics of the fluid velocity fields away from the probe in parallel and perpendicular directions for different values of $\beta _\circ$, $\nu$ and $r/a$. We observe that the fluid velocity does not undergo a slow relaxation process and the relaxation dynamics is dominated by shear relaxation times $\tau$. This difference can be explained by noting that the fluid, unlike the network, is incompressible, hence its dynamics is determined by shear hydrodynamic modes and not the compression modes. These results confirm further our observations thus far: i.e. for the network displacement fields the early/fast relaxation dynamics is controlled by $\tau$, and the slow relaxation is controlled by $\tau _D^r$, whereas the relaxation of fluid velocity is determined only by $\tau$. Another difference between VE and PVE models appears in the ratio of parallel to perpendicular displacements. For VE materials, $u_{VE}^\parallel /u_{VE}^\perp =M_{12}^{\parallel }/M_{12}^{\perp }=2$, whereas for a linear elastic material, the ratio is dependent on $\nu$, $u_{E}^\parallel /u_{E}^\perp ={4(1-\nu )}/{(3-4\nu )}$ (Levine & Lubensky Reference Levine and Lubensky2000), which changes from $2$ for $\nu =0.5$ to $\frac {8}{7}$ for $\nu =-1$. For PVE materials, the ratio changes over time from $2$ at early times to $u_{E}^\parallel /u_{E}^\perp$ at long times. Figure 3(a) shows these variations over time for $10\le r/a \le 50$ and $\beta _\circ =2$, while the inset shows the results for $2\le \beta _\circ \le 20$ and $r=50$. In both plots, $\nu =0.3$ and the time is made dimensionless by $\tau _D^r=r^2/D_q$. As can be seen, the results for different values of $r/a$ and $\beta _\circ$ all collapse to a single curve over the entire time domain.

Figure 3. (a) The ratio of parallel to perpendicular network displacements, $\varGamma _n=u^\parallel _n/u^\perp _n$, as a function of $t/\tau ^D_r$. Main plot: different values of distance $r$ from the probe, taking $\beta _\circ =2$ and $\nu =0.3$. Inset: different permeability $\beta _\circ$, taking $r/a=10$, $\nu =0.3$. (b) Variations of $\varTheta ={(\varGamma _n -\varGamma _{n}^{\infty })}/{(\varGamma _{n}^{\circ }-\varGamma _{n}^{\infty })}$ versus time for different Poisson ratios, taking $\beta _\circ =2$ and $r/a=10$. Inset: collection of all the data in the figure plotted as $\varTheta$ versus $t/\tau ^D_r$. All the data collapse to a single curve irrespective of the choice of $\nu$, $r$, $\beta _\circ$ and $\tau$.

Figure 4. Relaxation of displacements in the parallel (a) and perpendicular (b) directions to the applied force, for different values of permeability $\beta _\circ$, Poisson ratio $\nu$, and separation distances $r/a$. In all cases, the relaxation is dominated (down to $1-u^{(\parallel,\perp )}_f/u_f^\infty \approx 0.001$) by the shear relaxation time scale $\tau$.

In figure 3(b), we study the effect of $\nu$ on the variations of $\varTheta ={(\varGamma _n -\varGamma _{n}^\infty )}/{(\varGamma _{n}^{\circ }-\varGamma _{n}^\infty )}$ as a function of $t/\tau _D^r$ for different values of $\nu$, with $\beta _\circ =2$ and $r/a=10$, where $\varGamma _n (r,t)=u_n^{\parallel }/u_n^\perp$. The data collapse to a single curve. The inset of figure 3(b) shows the cumulative $\varTheta (t)$ versus $t/\tau _D^r$ data for different values of $r$, $\beta _\circ$ and $\nu$, collapsing to the same curve. All together, these results show that the relaxation dynamics of displacement (and mobility) ratios is determined entirely by the diffusion of network compressibility $\tau _D^r$. The experimental determination of this ratio in two-point MR provides another method to compute $D_q$ and thus $\beta _\circ$. In contrast, as we show in figure 5(a) of Appendix D, the ratio of fluid displacements in parallel and perpendicular directions remains nearly equal to $2$ for different choices of parameters.

Figure 5. (a) The ratio of parallel to perpendicular fluid displacements, $\varGamma _f=u^\parallel _f/u^\perp _f$, as a function of $t/\tau$. Main plot: different values of $r$, taking $\beta _\circ =2$ and $\nu =0.3$. Inset: different permeabilities, for $r/a=10$ and $\nu =0.3$. (b) The sphere relaxation $1-x(t)/x_\infty$ versus time for different values of viscosity, when $\nu =0.3$, $G=1$ and $\xi =1$. Inset plot represents the results as a function of rescaled time $t/\tau _D$, where $\tau _D={a^2 \xi (1-2\nu )}/{2G(1-\nu )}$ is the diffusion time scale of the compressibility field for distance $a$.

4. Summary

Increasing experimental evidence points to a PVE model as the most accurate description of the mechanics of the cell cytoskeleton (Mogilner & Manhart Reference Mogilner and Manhart2018). Biot's theory of poroelasticity has been used to describe these observations (Charras et al. Reference Charras, Yarrow, Horton, Mahadevan and Mitchison2005, Reference Charras, Coughlin, Mitchison and Mahadevan2008). Biot's theory neglects shear stresses, which are crucial for determining the mechanical response of cytoskeletal assemblies ($\phi \ll 1$) at short and intermediate time scales. To overcome this limitation, we include the shear stresses, which modifies the fluid momentum equation from Darcy to the Brinkman equation. We present the first general solutions of these modified linear PVE equations in spherical coordinates; the fluid phase is described as an incompressible linear VE fluid, and the network phase is modelled as a compressible linear VE material. Similar to Stokes flow microhydrodynamics (Kim & Karrila Reference Kim and Karrila2013) and linear elasticity (Gurtin Reference Gurtin1973), the linearity of the equations of motion allows development of a robust mathematical framework for describing the motion of inclusions moving in PVE materials, and the interior displacements generated in PVE materials by the motion of outer boundaries. These general solutions can be used to develop Faxén laws for a two-phase system, which relates the force (torque) on a sphere to the (angular) velocity fields of the fluid and network in general background flow. Faxén laws can then be used to develop numerical techniques similar to Stokesian dynamics (Fiore & Swan Reference Fiore and Swan2019), for modelling the motion of spherical inclusions in PVE materials, based on singularity solutions.

To demonstrate the utility of the general solutions, we studied the dynamics of a rigid sphere moving in a PVE material. Using the axisymmetric form of the general solutions, we derived closed-form solutions for the sphere's response function (hydrodynamic resistance), and fluid and network displacement fields. We showed that the displacement relaxation of the rigid sphere and the network displacements induced by it follow two distinct dynamics at short and long time scales. The dynamics at short time scales is dominated by hydrodynamic shear modes on the system, and thus can be captured using a VE model. The following slow relaxation dynamics could not be predicted using the VE model; instead, we showed that this dynamic is controlled by the diffusion time scale of network compressibility and the response of the fluid pressure to these local changes in network volume. These predictions are in agreement with the experimental observations, including the slow distance-dependent relaxation that follows the initial fast relaxation observed by atomic force microscopy experiments of Rosenbluth et al. (Reference Rosenbluth, Crow, Shaevitz and Fletcher2008). They are also in agreement with the results of Charras, Mitchison & Mahadevan (Reference Charras, Mitchison and Mahadevan2009) and Moeendarbary et al. (Reference Moeendarbary, Valon, Fritzsche, Harris, Moulding, Thrasher, Stride, Mahadevan and Charras2013), and the experimental and computational observations of slow propagation of hydrostatic pressure in blebs (Charras et al. Reference Charras, Yarrow, Horton, Mahadevan and Mitchison2005, Reference Charras, Coughlin, Mitchison and Mahadevan2008; Strychalski & Guy Reference Strychalski and Guy2016). Finally, we discussed how our results can be used directly to analyse the results of single-particle and two-point MR, and calculate fluid permeability in addition to the commonly measured VE properties of the system.

In studying the probe's motion, we assumed that the probe size is significantly larger than the mesh size of the network, which led to assuming no-slip BCs for both the fluid and network phases at the probe's surface ($r=a$). This is a reasonable assumption, since the mesh size for actin networks in the cell starts from 20 to 200 nm, and the typical probe size used in the experiments is at least 1000 nm, which is about at least $5$ times larger than the mesh size. On the other hand, it is now possible to include and trace genetically encoded multimeric nanoparticles of diameter 15–40 nm within the cell (Delarue et al. Reference Delarue2018). The particles that are smaller than the mesh size of the network can be used to study the effect of network topology and mechanics on the transport of small molecules and organelles through the fluid phase. The BCs for the network phase need to be modified when analysing the motion of particles when $\beta _\circ a\le 1$. Here, the more appropriate BCs for the network phase would be to assume that because the network and the probe are not in direct contact, no force is applied from the network to the probe, i.e. $\boldsymbol {\sigma }_n\boldsymbol {\cdot } \hat {\boldsymbol {n}}=\boldsymbol {0}$ at $r=a$ (Fu et al. Reference Fu, Shenoy and Powers2008; Diamant Reference Diamant2015). This BC can be incorporated easily within our general solution. Note that this modification in BCs changes qualitatively the time-dependent motion of the probe. In particular, the probe's response at long times will be dominated by VE properties of the fluid domain when $\beta _\circ a\le 1$, in contrast to the mechanics being determined by the network phase when $\beta _\circ a> 1$.