Oriented suspension mechanics with applications to flow linear dichrosim spectroscopy and pathogen detection

TL;DR

This paper develops a model for flow linear dichroism in suspension-based pathogen detection, showing how channel dimensions and flow oscillations influence particle alignment and signal sensitivity.

Contribution

It introduces a coupled Fokker-Planck and flow model for oriented suspension mechanics, optimizing parameters for enhanced linear dichroism signals in pathogen detection.

Findings

Increasing channel depth significantly improves the dichroism signal.

Oscillating flow can achieve similar alignment to steady flow at certain frequencies.

Optimal channel depth and pressure maximize particle alignment.

Abstract

Flow linear dichroism is a biophysical spectroscopic technique that exploits the shear-induced alignment of elongated particles in suspension. Motivated by the broad aim of optimising the sensitivity of this technique, and more specifically by a handheld synthetic biotechnology prototype for waterborne-pathogen detection, a model of steady and oscillating pressure-driven channel flow and orientation dynamics of a suspension of slender microscopic fibres is developed. The model couples the Fokker-Planck equation for Brownian suspensions with the narrow channel flow equations, the latter modified to incorporate mechanical anisotropy induced by the particles. The linear dichroism signal is estimated through integrating the perpendicular components of the distribution function via an appropriate formula which takes the bi-axial nature of the orientation into account. For the specific application of pathogen detection via binding of M13 bacteriophage, it is found that increases in the channel depth are more significant in improving the linear dichroism signal than increases in the channel width. Increasing channel depth to $2~$mm and pressure gradient to $5 \times 10^4~$Pa/m essentially maximises the alignment. Oscillating flow can produce nearly equal alignment to steady flow at appropriate frequencies, which has significant potential practical value in the analysis of small sample volumes.

Figures

Tables

Table 1. Table 1: Summary of the standard parameter set split into three categories: (a) the parameters that vary throughout the results, (b) the fixed parameters and physical constants and (c) the dimensionless groups.

Parameter	Description	Standard value	Range	Units
$h^{\ast }$	Channel half-depth	$3\times {10}^{-4}$	$2\times {10}^{-4}$ – $4\times {10}^{-3}$	m
$W^{\ast }$	Channel half-width	$1\times {10}^{-3}$	$0.4\times {10}^{-3}$ – $2.5\times {10}^{-3}$	m
$n_d^{\ast }$	Number density	$7.33\times {10}^{17}$	$1\times {10}^{16}$ – $1\times {10}^{18}$	phage/m3
$G^{\ast }$	Pressure gradient	$4.4\times {10}^3$	$0.4\times {10}^3$ – $10\times {10}^3$	Pa/m
${\omega }^{\ast }$	Frequency of oscillations	$10$	$5$ – $70$	s-1

Table 2. Table 2: The spatially averaged orientation parameter for oscillating flow and a range of frequency of oscillations, pressure gradient and channel depth. Throughout these results, n d ∗ = 1 × 10 18 superscript subscript 𝑛 𝑑 1 superscript 10 18 n_{d}^{*}=1\times 10^{18}~{} phage/m 3 and W ∗ = 2 × 10 − 3 superscript 𝑊 2 superscript 10 3 W^{*}=2\times 10^{-3}~{} m.

$G^{\ast }$	$h^{\ast }$	$\overline{S}$	$\overline{S}$
		(${\omega }^{\ast }=10$s-1)	(${\omega }^{\ast }=100$s-1)
$2\times {10}^3$	$2\times {10}^{-4}$	$0.130$	$0.042$
	$4\times {10}^{-4}$	$0.193$	$0.031$
$4\times {10}^3$	$2\times {10}^{-4}$	$0.241$	$0.110$
	$4\times {10}^{-4}$	$0.316$	$0.081$
$6\times {10}^3$	$2\times {10}^{-4}$	$0.309$	$0.165$
	$4\times {10}^{-4}$	$0.381$	$0.122$
$8\times {10}^3$	$2\times {10}^{-4}$	$0.353$	$0.207$
	$4\times {10}^{-4}$	$0.424$	$0.154$

Table 3. Table 3: The spatially averaged orientation parameter for oscillating flow and a range of channel width, number density and frequencies of oscillations. Throughout these results, G ∗ = 8 × 10 3 superscript 𝐺 8 superscript 10 3 G^{*}=8\times 10^{3}~{} Pa/m and n d ∗ = 1 × 10 18 superscript subscript 𝑛 𝑑 1 superscript 10 18 n_{d}^{*}=1\times 10^{18}~{} phage/m 3 .

$W^{\ast }$	$h^{\ast }$	$\overline{S}$	$\overline{S}$
		(${\omega }^{\ast }=10$s-1)	(${\omega }^{\ast }=100$s-1)
$0.5\times {10}^{-3}$	$2\times {10}^{-4}$	$0.302$	$0.227$
	$4\times {10}^{-4}$	$0.406$	$0.217$
$2\times {10}^{-3}$	$2\times {10}^{-4}$	$0.287$	$0.212$
	$4\times {10}^{-4}$	$0.373$	$0.173$

Equations

$${\rho }^{\ast }\left(\frac{\partial {\mathbf{u}}^{\ast }}{\partial t^{\ast }}+({\mathbf{u}}^{\ast }\cdot {\nabla }^{\ast }){\mathbf{u}}^{\ast }\right)$$

$${\nabla }^{\ast }\cdot {\mathbf{u}}^{\ast }$$

$${\mathbf{u}}^{\ast }(x^{\ast },\pm W^{\ast },z^{\ast },t^{\ast })=0,{\mathbf{u}}^{\ast }(x^{\ast },y^{\ast },\pm h^{\ast },t^{\ast })=0.$$

$$\frac{\partial \psi }{\partial t^{\ast }}+{\nabla }_{\mathbf{p}}\cdot ({\mathbf{\Omega }}^{\ast }\psi )=D_r^{\ast }{\nabla }_{\mathbf{p}}^2\psi ,$$

$$\mathbf{p}=(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta ),$$

$${\mathbf{\Omega }}^{\ast }=(\text{I}-\mathbf{pp})\cdot ({\alpha }_0{\mathbf{e}}^{\ast }+{\mathbf{\omega }}^{\ast })\cdot \mathbf{p},$$

$$\psi (\varphi ,\theta )$$

$${\int }_s\psi d\mathbf{p}=1,$$

$${\mathbf{\sigma }}^{\ast }={\mathbf{\sigma }}_I^{\ast }+{\mathbf{\sigma }}_D^{\ast }+{\mathbf{\sigma }}_P^{\ast },$$

$${\mathbf{\sigma }}_I^{\ast }=-P^{\ast }\text{I}+2{\mu }^{\ast }{\mathbf{e}}^{\ast },$$

$${\mathbf{\sigma }}_D^{\ast }=$$

$$\displaystyle\begin{split}\boldsymbol{\sigma}_{P}^{*}=&4\,\mu^{*}\,\Phi\Bigg{[}\alpha_{2}\,{\boldsymbol{e}}^{*}\mathrel{\mathop{\mathchar 58\relax}}\int_{s}\boldsymbol{p\,p\,p\,p}\,\psi\operatorname{d\!}{\bf p}+\alpha_{3}\left({\boldsymbol{e}}^{*}\cdot\int_{s}\boldsymbol{p\,p}\,\psi\operatorname{d\!}{\boldsymbol{p}}+\int_{s}\boldsymbol{p\,p}\,\psi\operatorname{d\!}{\boldsymbol{p}}\cdot{\boldsymbol{e}}^{*}\right)\\ &\hskip 85.35826pt+\alpha_{4}\,{\boldsymbol{e}}^{*}\int_{s}\psi\operatorname{d\!}{\boldsymbol{p}}+\alpha_{5}\,{\boldsymbol{I}}\,{\boldsymbol{e}}^{*}\mathrel{\mathop{\mathchar 58\relax}}\int_{s}\boldsymbol{p\,p}\,\psi\operatorname{d\!}{\boldsymbol{p}}\Bigg{]},\end{split}$$

$$\begin{array}{r}{u}^{\ast }=\frac{G^{\ast }{h^{\ast }}^2}{{\mu }^{\ast }}u,v^{\ast }=\frac{\delta G^{\ast }{h^{\ast }}^2}{{\mu }^{\ast }}v,w^{\ast }=\frac{\delta G^{\ast }{h^{\ast }}^2}{{\mu }^{\ast }}w,\\ x^{\ast }=\frac{h^{\ast }}{\delta }x,y^{\ast }=h^{\ast }y,z^{\ast }=h^{\ast }z,P^{\ast }=\frac{G^{\ast }{h}^{\ast }}{\delta }P,\end{array}$$

$$\begin{array}{rl}S=& ⟨p_{x^{\ast }}^2⟩-⟨p_{y^{\ast }}^2⟩,\\ =& ⟨{\sin }^2\theta {\cos }^2\varphi ⟩-⟨{\sin }^2\theta {\sin }^2\varphi ⟩,\end{array}$$

$$⟨\bullet ⟩={\int }_0^{2\pi }{\int }_0^{\pi }\bullet \psi \sin \theta d\theta d\varphi .$$

$$S=\frac{1}{2}\left(3{\int }_s{\cos }^2{\theta }_S\psi d\mathbf{p}-1\right),$$

$${\theta }_S=\arccos (|\sin \theta \cos \varphi |).$$

$$\frac{1}{6}\Lambda (\psi )=P_e{\rm Y}(\psi ),$$

$$\Lambda (\psi )$$

$${\rm Y}(\psi )$$

$$\psi ({\varphi }^{\prime },{\theta }^{\prime })=\sum _{n=0}^N\sum _{m=0}^n\left(A_{0n}^m{P}_n^m(\cos {\theta }^{\prime })\cos (m{\varphi }^{\prime })+A_{1n}^m{P}_n^m(\cos {\theta }^{\prime })\sin (m{\varphi }^{\prime })\right),$$

$$\Lambda (P_n^m(\cos {\theta }^{\prime })\cos (m{\varphi }^{\prime }))$$

$$\Lambda (P_n^m(\cos {\theta }^{\prime })\sin (m{\varphi }^{\prime }))$$

$${\rm Y}(P_n^m(\cos {\theta }^{\prime })\cos (m{\varphi }^{\prime }))$$

$${\rm Y}(P_n^m(\cos {\theta }^{\prime })\sin (m{\varphi }^{\prime }))$$

$$\frac{q(q+1)}{6}{A}_{0q}^p$$

$$\frac{q(q+1)}{6}{A}_{1q}^p$$

$$\begin{array}{rl}-1=& \frac{\partial }{\partial y}(\left\{1+4\Phi \left[{\alpha }_2{\int }_s{p}_1^2{p}_2^2\psi d\mathbf{p}+\frac{{\alpha }_3}{2}{\int }_s(p_1^2+p_2^2)\psi d\mathbf{p}+\frac{{\alpha }_4}{2}\right]\right\}\frac{\partial u}{\partial y}\\ & +4\Phi \left[{\alpha }_2{\int }_s{p}_1^2{p}_2{p}_3\psi d\mathbf{p}+\frac{{\alpha }_3}{2}{\int }_s{p}_2{p}_3\psi d\mathbf{p}\right]\frac{\partial u}{\partial z})\\ +& \frac{\partial }{\partial z}(\left\{1+4\Phi \left[{\alpha }_2{\int }_s{p}_1^2{p}_3^2\psi d\mathbf{p}+\frac{{\alpha }_3}{2}{\int }_s(p_1^2+p_3^2)\psi d\mathbf{p}+\frac{{\alpha }_4}{2}\right]\right\}\frac{\partial u}{\partial z}\\ & +4\Phi \left[{\alpha }_2{\int }_s{p}_1^2{p}_2{p}_3\psi d\mathbf{p}+\frac{{\alpha }_3}{2}{\int }_s{p}_2{p}_3\psi d\mathbf{p}\right]\frac{\partial u}{\partial y})\\ +& \frac{2\Phi {\alpha }_r}{P_G}\left(\frac{\partial }{\partial y}{\int }_s{p}_1{p}_2\psi d\mathbf{p}+\frac{\partial }{\partial z}{\int }_s{p}_1{p}_3\psi d\mathbf{p}\right),\end{array}$$

$$u(\pm \overline{W},z)=0,u(y,\pm 1)=0,$$

$$\frac{\partial P^{\ast }}{\partial x^{\ast }}=G^{\ast }\exp (i{\omega }^{\ast }{t}^{\ast }),$$

$$\frac{\partial \psi }{\partial \tau }=\frac{1}{6}\Lambda (\psi )-P_e(\tau ){\rm Y}(\psi ),$$

$$\psi ({\varphi }^{\prime },{\theta }^{\prime },t)=\sum _{n=0}^N\sum _{m=0}^n\left(A_{0n}^m(t)P_n^m(\cos {\theta }^{\prime })\cos (m{\varphi }^{\prime })+A_{1n}^m(t)P_n^m(\cos {\theta }^{\prime })\sin (m{\varphi }^{\prime })\right).$$

$$\frac{d{A}_{0q}^p}{dt}$$

$$\frac{d{A}_{1q}^p}{dt}$$

$$\begin{array}{rl}{\alpha }^2\frac{\partial u}{\partial t}=& -\exp (it)+\frac{2\Phi {\alpha }_r}{P_G}\left(\frac{\partial }{\partial y}{\int }_s{p}_1{p}_2\psi d\mathbf{p}+\frac{\partial }{\partial z}{\int }_s{p}_1{p}_3\psi d\mathbf{p}\right)\\ +& \frac{\partial }{\partial y}(\left\{1+4\Phi \left[{\alpha }_2{\int }_s{p}_1^2{p}_2^2\psi d\mathbf{p}+\frac{{\alpha }_3}{2}{\int }_s(p_1^2+p_2^2)\psi d\mathbf{p}+\frac{{\alpha }_4}{2}\right]\right\}\frac{\partial u}{\partial y}\\ & +4\Phi \left[{\alpha }_2{\int }_s{p}_1^2{p}_2{p}_3\psi d\mathbf{p}+\frac{{\alpha }_3}{2}{\int }_s{p}_2{p}_3\psi d\mathbf{p}\right]\frac{\partial u}{\partial z})\\ +& \frac{\partial }{\partial z}(\left\{1+4\Phi \left[{\alpha }_2{\int }_s{p}_1^2{p}_3^2\psi d\mathbf{p}+\frac{{\alpha }_3}{2}{\int }_s(p_1^2+p_3^2)\psi d\mathbf{p}+\frac{{\alpha }_4}{2}\right]\right\}\frac{\partial u}{\partial z}\\ & +4\Phi \left[{\alpha }_2{\int }_s{p}_1^2{p}_2{p}_3\psi d\mathbf{p}+\frac{{\alpha }_3}{2}{\int }_s{p}_2{p}_3\psi d\mathbf{p}\right]\frac{\partial u}{\partial y}),\end{array}$$

Full text

1 Introduction

Suspensions of particles in liquid or gas are found throughout the natural world and in many industrial processes; including blood, particulate-laden air, algae in open water or in bioconvection experiments, semen samples, and industrial applications such as emulsions in food and cosmetics – these examples and more are elaborated in references [32, 21, 7, 20]. We will consider flowing suspensions of microscopic elongated rod-like particles, as a means to understand and improve flow linear dichroism spectroscopy, with particular focus on a specific technological application. The core idea of flow linear dichroism spectroscopy is that elongated particles undergo shear-induced rotation, which has the effect of concentrating their orientation in the direction of flow. The difference in absorbance of light polarised parallel, and perpendicular to, the flow direction can be used to reveal structural information; further information and wider context is given in references [39, 10, 6, 47]. The application of mathematical modelling to understand and quantify molecular orientation in flow linear dichroism was addressed by McLachlan et al. [40], building on classical oriented suspension mechanics [3, 49]. This manuscript will generalise this work to the significantly more complex problem of a non-homogeneous shear environment of pressure-driven, and potentially time-varying, channel flow, and will moreover focus on a specific technological application.

The specific application is a prototype handheld device, developed by Linear Diagnostics Ltd (LD) designed to detect waterborne pathogens in fluids. The analyte is mixed with a reagent containing a synthetic biology micron-length fibre based on M13 bacteriophage. In the absence of the pathogen, the fibres align in the flow field and therefore exhibit linear dichroism. The M13 are engineered with antibodies enabling binding to specific pathogens, therefore in the presence of the pathogen the alignment is disrupted, resulting in a reduced linear dichroism signal. The technique is described further in reference [43]. To improve the sensitivity of the device it is desirable to optimise the signal to noise ratio, which can be achieved by increasing the overall alignment in the system. To analyse the system theoretically involves solving a coupled pressure-driven steady or oscillating flow and alignment problem. While we present results for a specific biotechnology application, the computational framework is generic to any channel flow of dilute suspensions.

The coupling between the flow and particle orientation is important to understand the dynamics of these suspensions. Jeffrey [29] carried out a landmark study of the dynamics of suspended elongated particles in viscous flow, in particular shear induced alignment. Moreover, following seminal theoretical work [15, 35, 2, 45], it has been established how the aspect ratio of particles affects emergent rheology. In the related context of suspensions of actively-swimming oriented particles, Pedley & Kessler [44] described a coupled model of orientation distribution, governed by a Fokker-Planck equation, and fluid flow, given by modifications to the Navier-Stokes equations proportional to the ensemble averages of orientation. These early studies form the theoretical framework for our model of homogeneous suspensions of semi-rigid M13 bacteriophage via the approximation of elongated, axisymmetric, rigid Brownian rods in the limit of infinite aspect ratio. Other relevant theoretical studies include work on homogeneous shear [49, 30, 34], time-dependent shear [22], turbulent channel flow [37, 38, 41] and pressure-driven channel flow [16]. More recently, Taylor-Couette [26] and Rayleigh-Bénard [27] stability of perfectly aligned suspensions has been theoretically modelled via the transversely isotropic fluid of Ericksen [15, 19]. Similar models have also been used to investigate fibrous biological systems such as growth in plant cell walls [14], propulsion in aligned cervical mucus [8, 9] and the extracellular matrix [13].

This paper solves and compares two 3D pressure-driven flow problems for an oriented suspension in a narrow rectangular channel geometry: steady flow, which models the continuously pumped loop used in existing spectroscopic devices, including the LD technology, and a potential novel oscillatory system, in which a smaller volume of analyte is oscillated back and forth. For practical particle concentrations we have $n_{d}^{*}{a^{*}}^{3}<1$, where $n_{d}^{*}$ is particle number density and $a^{*}$ is the particle length, and hence the suspension is at the upper end of the dilute range [11]. The channel dimensions, pressure gradient, particle number density and frequency of oscillations are investigated as factors to improve alignment, and thus signal, in both flow types and to determine the viability of an oscillatory system for aligning particles.

The coupled orientation and flow model will be solved numerically by iterative coupling. Mathematical modelling of these suspensions is computationally challenging due to the additional independent variables associated with the particle degrees of freedom and the coupling between particle dynamics, velocity gradients and rheology. Rational simplification of the flow problem via lubrication theory, the application of a spectral method [49] and multicore parallelisation of the array of spatially-local orientation problems, will be shown to enable practical solution with workstation hardware.

The manuscript is organised as follows: the governing equations for the system, including the Navier-Stokes and Fokker-Planck equations, are summarised in section 2. The steady flow model is presented in section 3 and the oscillatory problem in section 4. Results for both flow problems are presented in section 5 and discussed and compared in section 6.

2 Summary of equations governing dilute suspensions of elongated particles

Consider a 3D rectangular channel of width $2W^{*}$, depth $2h^{*}$ and length scale $L^{*}$, where $h^{*},\,W^{*}\ll L^{*}$. The axis origin is at the centre of the channel and it is assumed that the flow direction induced by pressure gradient $G^{*}$, or the molecular orientation axis, is the $x^{*}$-direction (figure 1 (a)). The behaviour of dilute suspensions is approximated by altering the constitutive equation for the stress tensor $\boldsymbol{\sigma}^{*}$ [23, 24]; the motivating problem is at the upper end of the dilute regime (see section 5). This constitutive equation, along with the coupled Fokker-Planck equation governing fibre behaviour [49], are briefly detailed now.

The dimensional Navier-Stokes and incompressibility equations are

[TABLE]

where $\rho^{*}$ is density, $\boldsymbol{u}^{*}(\boldsymbol{x}^{*},t^{*})$ is the velocity vector for spatial coordinate $\boldsymbol{x}^{*}$ and time $t^{*}$ and $\nabla^{*}$ is differentiation in Cartesian space. The asterisk notation represents dimensional variables. The boundary conditions, are given by no-slip and no-penetration at the walls,

[TABLE]

The boundary conditions for $x^{*}$ have not been defined here as they are specified for each flow system.

To determine the bulk stresses in the fluid, consider a homogeneous suspension of particles, dilute enough to ensure particle-particle interactions are negligible [49, 40]. Further assume the suspension is spatially homogeneous, then the suspension configuration can be described by the orientation distribution function $\psi(\boldsymbol{x}^{*},\boldsymbol{p},t^{*})$, satisfying the Fokker-Planck equation,

[TABLE]

where $\boldsymbol{p}$ is the unit orientation vector,

[TABLE]

for $\theta\in[0,\pi]$ and $\phi\in[0,2\pi)$ shown in figure 1 (b), $\nabla_{\boldsymbol{p}}$ represents differentiation in orientation space and $D_{r}^{*}$ is the rotational diffusion constant. The evolution of the director $\boldsymbol{p}$ is given by

[TABLE]

where $\boldsymbol{e}^{*}=(\nabla^{*}\boldsymbol{u}^{*}+\nabla^{*}{\boldsymbol{u}^{*}}^{T})/2$ is the rate of strain tensor, $\boldsymbol{\omega}^{*}=(\nabla^{*}\boldsymbol{u}^{*}-\nabla^{*}{\boldsymbol{u}^{*}}^{T})/2$ is the vorticity tensor and $\alpha_{0}=(r^{2}-1)/(r^{2}+1)$, for particle aspect ratio $r$, is a measure of slenderness [29]. The orientation distribution function must satisfy periodicity and the normalisation condition,

[TABLE]

where $s$ is the unit surface of a sphere.

Following Pedley and Kessler [44] (and the prior work of Hinch and Leal [24]), the enhanced stress tensor for a suspension of inactive particles is

[TABLE]

where the Newtonian contribution is

[TABLE]

for pressure $P^{*}$, identity matrix I and viscosity $\mu^{*}$. Here, $\boldsymbol{\sigma}_{D}^{*}$ is the extra stress due to particle rotation and $\boldsymbol{\sigma}_{P}^{*}$ represents the interaction of the particles with the fluid, given by

[TABLE]

and where $\Phi$ is the volume fraction of particles. The constants $\alpha_{i}$ ($i=2\dots 5,\,r$) relate to the aspect ratio of the particle and can be represented in terms of ellipsoidal integrals [2, 29], which are detailed in equations (S1)–(S5) in supplemental materials section S1. For reference, Holloway et al [25] have summarised this system of equations for particle suspensions in concise notation.

Since the width and depth of the channel are comparable and much smaller than the length scale, we use a lubrication approximation, introducing the dimensionless parameter $\delta=h^{*}/L^{*}\ll 1$ to allow detailed analysis of the significant hydrodynamic interactions. The variables are scaled via

[TABLE]

where $G^{*}$ is a reference pressure gradient. For steady flow, $G^{*}$ is a complete pressure gradient and for oscillating flow it is determined from the amplitude of the oscillating pressure gradient. The full dimensionless models for steady and oscillatory flow are described in the relevant sections.

Finally, the degree of orientation in a system can be described by an orientation parameter, $S$, calculated as

[TABLE]

for a 3D suspension where

[TABLE]

Note that this is different from the definition used by [40],

[TABLE]

where

[TABLE]

The expression (16) relates only to uni-axial samples such as planar membranes in which $\psi=\psi(\theta_{S},t^{*})$ [17, 47], rather than flow spectroscopy in which the orientation distribution is bi-axial, i.e. $\psi=\psi(\theta,\phi,t^{*})$. McLachlan’s[40] use of this formula is therefore inaccurate in this context.

Equations (1)–(12) constitute the full dimensional model for variables $u^{*}(\boldsymbol{x}^{*},t^{*})$ and $\psi(\phi,\theta,t^{*})$ which will be solved numerically. The steady flow system is described in section 3 and the oscillatory system in section 4.

3 Steady flow

The first problem we consider is flow forced by a constant pressure gradient $G^{*}$. Section 3 3.1 contains the Fokker-Planck equation and section 3 3.2 contains the model for the flow. The full coupled problem is solved by iterating between the numerical methods described in section 3 3.1 and section 3 3.2.

3.1 Model for the orientation distribution function

The particles in the flow are small enough that each particle at a point in space is subject to an approximately linear velocity field in the local frame, and hence a constant shear rate $\dot{\gamma}^{*}=\sqrt{(\nabla^{*}\boldsymbol{u}^{*})^{2}+(\nabla^{*}{\boldsymbol{u}^{*}}^{T})^{2}}$. Therefore the orientation distribution of the bacteriophage is formulated at each point in space independently; fluid flow influences orientation distribution only via the local Péclet number. We describe the model and its numerical discretisation below. Following Strand, Kim and Karrila [49] (see also Bird et al [5]), assuming $\alpha_{0}=1$ and investigating steady state behaviour, the Fokker-Planck equation reduces to

[TABLE]

where $P_{e}=\dot{\gamma}^{*}/6D_{r}^{*}$ is the local Péclet number and

[TABLE]

where $\theta^{\prime}$ and $\phi^{\prime}$ describe the local coordinates at each point in space, which are calculated by rotating $\phi$ and $\theta$ away from the lab frame by an angle $\beta=\text{acos}(|\partial u^{*}/\partial y^{*}|/\dot{\gamma}^{*})$.

The Fokker-Planck equation (17) is spatially discretised via spherical harmonics [5, 40, 49]. Assume the solution takes the form

[TABLE]

where $P_{n}^{m}$ are the associated Legendre polynomials. By substituting (20) into equation (17), $\Lambda$ and $\Upsilon$ are given by

[TABLE]

The first two results follow from the definition of a spherical harmonic, i.e. the solution to the spherical part of Laplace’s equation on a sphere; the other two are a linear combination of spherical harmonics where the seven non-zero constants $a_{n,q}^{m,p}$, reproduced from Bird [5] and Strand, Kim & Karrila [49], are given in supplemental materials section S2. The resulting system of equations for $A_{0\,n}^{m}$ and $A_{1\,n}^{m}$, found by equating coefficients of sine and cosine, is

[TABLE]

for integer $N$, $q=0,2,\ldots,N$ and $p=0,2,\ldots,q$. Here, $A_{0\,0}^{0}=1/4\pi$ to satisfy the normalisation condition, and we have $A_{1\,n}^{0}=0$ for all $n$. Due to particle symmetry, if $n$ or $m$ are odd then $A_{0\,n}^{m}=A_{1\,n}^{m}=0$.

The system (25)–(26) is solved numerically by setting up a matrix of equations for $A_{0\,q}^{p}$ and $A_{1\,q}^{p}$ and using a direct solver (backslash) in Matlab.

3.2 Model for the velocity

Under the lubrication theory approximation, the flow problem can be expressed entirely as an elliptic partial differential equation for $u(y,z)$ where the coefficients depend on moments of the orientation distribution $\psi(\phi,\theta,y^{*},z^{*})$.

Scaling the system, via (13), neglecting terms $\mathcal{O}(\delta)$ and smaller, the continuity equation (2) is unchanged and the Navier-Stokes equations (1) reduce to

[TABLE]

for global Péclet number $P_{G}=G^{*}h^{*}/D_{r}^{*}\mu^{*}$. The $y$– and $z$–components of the dimensionless Navier-Stokes equations enforce $P=P(x)$ and so the pressure gradient reduces to $-1$. The boundary conditions (3) become

[TABLE]

where $\overline{W}=W^{*}/h^{*}$. Equation (27), along with boundary conditions (28), are solved numerically using a finite difference scheme and constructing a matrix of equations from the resulting problem. The numerical solution to the Newtonian and suspension models are described in supplemental materials section S3 and the iterative process is described in section S4.

To optimise the efficiency of the above processes, the convergence of the spherical harmonic solution and the finite difference approximation are investigated for a range of spherical harmonic modes $N$ and finite difference grid spacings $Y$ and $Z$. The full details and related figures are presented in supplemental materials section S3; we find that $N=10$, $Y=150$ and $Z=120$ are sufficient for convergence. We will discuss the results in section 5.

4 Oscillatory flow

The next flow system of interest consists of an oscillating pressure gradient where the length is still much larger than its width and depth but a smaller volume of fluid oscillates back and forth. The pressure gradient is now defined as

[TABLE]

for frequency of oscillations $\omega^{*}$. The dimensional scalings are given by equation (13), and the time is scaled against the frequency of oscillations as $t^{*}=t/\omega^{*}$.

The model for the time dependent Fokker-Planck equation is given in section 4 4.1 and the model for the velocity profile is given in section 4 4.2. The full coupled problem comes from iterating numerically between 4 4.1 and 4 4.2 over each time step.

4.1 Model for the orientation distribution

Assuming the particles are small enough that at each point in the channel they are subjected to a spatially linear velocity field and thus a shear rate dependent only on time $\dot{\gamma}^{*}=\dot{\gamma}^{*}(t^{*})$; the orientation distribution is again formulated at each point in space independently. Following section 3 3.1, the Fokker-Planck equation for shear rate $\dot{\gamma}(t)$ is

[TABLE]

where $P_{e}=\dot{\gamma}^{*}/6D_{r}^{*}$ is the local Péclet number, $t^{*}=\tau/6D_{r}^{*}$ and $\Lambda$ and $\Upsilon$ are given by (18) and (19) respectively.

The orientation distribution is discretised spatially in the basis of spherical harmonics,

[TABLE]

This is then substituted into (30), where $\Lambda$ and $\Upsilon$ are given by (21)-(24). Equating coefficients of sine and cosine, the resulting system of ordinary differential equations for $A_{0\,n}^{m}$ and $A_{1\,n}^{m}$ is

[TABLE]

for integer $N$, $q=0,2,\ldots,N$ and $p=0,2,\ldots,q$. Here, $A_{0\,0}^{0}(t)=1/4\pi$ to satisfy the normalisation condition, $A_{1\,n}^{0}=0$ for all $n$ and $A_{0\,n}^{m}=A_{1\,n}^{m}=0$ for odd $n,\,m$.

The ordinary differential equations (32) and (33) are solved numerically using a second order explicit Improved Euler method by first writing the equations in matrix form as described in section 3 3.1.

4.2 Model for the velocity

Using a lubrication theory approximation, the velocity problem can be expressed as a time dependent partial differential equation where coefficients of $u(y,z,t)$ depend upon moments of the orientation distribution function $\psi(\phi,\theta,y^{*},z^{*},t^{*})$.

Scaling the Navier-Stokes and continuity equations, (1) and (2) respectively, retaining terms larger than $\mathcal{O}(\delta)$, the continuity equation is unchanged and the dimensionless Navier-Stokes equations become

[TABLE]

for global Péclet number $P_{G}$ and Womersley number $\alpha^{2}=\omega^{*}{h^{*}}^{2}\rho^{*}/\mu^{*}$. It is assumed the system is initially at rest, accounted for by altering the dimensionless pressure gradient to $\partial P/\partial x=\exp(i(t+\pi/2))$.

Equation (34) is solved via an Alternating Direction Implicit (ADI) method in which the time step is split in half and one spatial variable is discretised in each time step [12]; mixed derivatives are treated explicitly throughout this analysis. As in section 3, a matrix equation is constructed for each time step and solved numerically; the details of this method can be found in supplemental materials section S5 and a description of the iterative procedure is described in S6.

A numerical convergence study was undertaken, detailed in supplemental materials section S5. As in section 3, $Y=150$ and $Z=120$ are satisfactory in the finite difference scheme and $T=1200$ time steps ensure convergence.

5 Results

The velocity, orientation parameter and linear dichroism signal are compared for a range of particle number density, channel dimensions and pressure gradient. We will refer to the standard parameter set for these quantities as: channel half-depth $h^{*}=3\times 10^{-4}~{}$m, half-width $W^{*}=1\times 10^{-3}~{}$m, fixed pressure gradient $G^{*}\approx 4.4\times 10^{3}~{}$Pa/m, particle number density $n_{d}^{*}=7.33\times 10^{17}~{}$phage/m3 and frequency of oscillations $\omega^{*}=10~{}$s*-1* (table 1 (a)). Using a particle length of $a^{*}=800~{}$nm, the value $n_{d}^{*}{a^{*}}^{3}=0.3753<1$ and so we argue that the approximation of a dilute regime is reasonable. The dimensionless groups are detailed in table 1 (c) and are subject to change depending on their composition of dimensional values; we do not directly vary the dimensionless groups. All other parameters are unchanged (table 1 (b)).

The notation $\overline{\,\bullet\,}$ and $\overline{\,\bullet\,}^{y^{*}}$ represent the spatial average and the width average for steady flow respectively,

[TABLE]

The linear dichroism signal is assumed proportional to the particle number density and is calculated as the root mean square of the orientation parameter across the depth of the channel. The root mean square of the orientation parameter will be denoted as $S_{rms}$ and the signal is given by

[TABLE]

where $N_{A}^{*}$ is the Avogadro constant ($\approx 6\times 10^{23}$ phage/mol), $m_{w}^{*}$ is the molecular weight of M13 bacteriophage ($m_{w}^{*}\approx 2\times 10^{7}~{}$g/mol from Linear Diagnostics Ltd) and $\epsilon^{*}$ is the extinction coefficient (m2/mol), which is a measure of the absorption of light through the sample [10]. The extinction coefficient is taken as $0.38~{}$m2/g at light wavelength $269~{}$nm [1, 18, 42]. A number of the parameters in table 1 and their calculation are discussed in supplemental materials section S7. As a direct comparison with steady flow, the maximum value of the orientation parameter in the oscillating channel is presented.

The computational walltime for this problem was large, especially for oscillating flow. The numerical solution was parallelised using ‘parfor’ in Matlab by splitting the solution to the Fokker-Planck equation, solved at each point in space, across five cores on a workstation (2017 Lenovo Thinkstation P710; Intel(R) Xeon(TM) E5-2646 CPU @ 2.40 GHz; 128 GB 2400 MHz RDIMM RAM). The runtime varied between $14$ and $20$ hours for a single parameter tuple; $175$ tuples were calculated for the oscillating results.

5.1 Steady pressure gradient

Figure 2 depicts the velocity profile, orientation parameter and linear dichroism signal for steady flow through the cross section of the channel, with the standard parameter set. As may be expected from the Newtonian solution[46], the velocity $u^{*}$ increases from zero at the channel walls to a maximum at the centre point (figure 2 (a)), corresponding to a minimum in the orientation parameter $S$ (figure 2 (b)). Associated with the large shear rate, the linear dichroism signal is largest near the channel walls and decreases to a local minimum at the centre (figure 2 (c)); there is an increase at the middle width of the channel ($y^{*}=0$) where larger shear near the walls at $z^{*}=\pm h^{*}$ increases the orientation parameter. Local minima also exist at the channel corners.

The orientation distribution function $\psi$ is displayed for varying points in the channel cross section (figure 3). At the centre of the channel where the shear rate is zero (figure 3 (b)) there are no biasing effects on the particles and so $\psi$ is small. Close to the walls at $y^{*}=\pm W^{*}$m the behaviour of $\psi$ is similar (figures 3 (c) and (d)); the particles prefer to align in the flow direction, for orientation angles $\phi=0,\,\pi$ and $2\pi$ and $\theta=\pi/2$. Close to the corner of the channel this behaviour is more prominent. Near $z^{*}=h^{*}$m (figure 3 (e)), the M13 bacteriophage are more likely to be aligned at $\phi=0$ and $\pi$ (note that the particles have no polarity so the two angles are equivalent) and $\theta=\pi/2$.

The flow profile accounting for a particle suspension is compared with the Newtonian profile $u_{N}^{*}$ for the same flow setup in figure 4. The actual difference shows that the profile for both flow systems is qualitatively the same; the suspension model produces a lower velocity however the change is small (figure 4 (a)). To confirm that this change is small, the reduction in velocity relative to the magnitude of the Newtonian velocity is also calculated in figure 4 (b) and the change is on the order of $10^{-3}$.

We compare the width averaged orientation parameter calculated in this work (solid line) with equation (16) (dashed line); while both definitions show maxima and minima of the orientation parameter at the same point in the channel, McLachlan’s definition for the orientation parameter [40] over-predicts the orientation in general.

In figure 6, the impact of the pressure gradient $G^{*}$ and particle number density $n_{d}^{*}$ on the spatially averaged orientation parameter $\overline{S}$ and linear dichroism signal $\overline{LD}$ is examined. Increasing the number density reduces the orientation parameter less than $0.5\%$ over two orders of magnitude (results not shown here) and increasing the pressure gradient produces a (non-linear) increase in $\overline{S}$, seen in figure 6 (a). However, increasing the number density significantly increases the linear dichroism signal due to its linear relationship with number density (figure 6 (b)).

The spatially averaged orientation parameter is calculated for a range of channel width and depth in figure 7; the other parameters, including pressure gradient, remain fixed. The channel depth has a large impact on $\overline{S}$; for a deeper channel, the orientation parameter is larger. The channel width has limited effect, unless the two dimensions are comparable (figure 7 (b)). For channel widths below $W^{*}=1\times 10^{-3}~{}$m increasing the channel width increases the orientation parameter; this effect reverses as $W^{*}$ grows. Finally, the upper asymptote of the spatially averaged orientation parameter is investigated for large pressure gradient and channel depth (figure 8). In both cases, $\overline{S}$ quickly tends towards a constant value, $\overline{S}\approx 0.74$.

5.2 Oscillating pressure gradient

The velocity, orientation parameter and linear dichroism signal for oscillatory flow with the standard parameter set are shown in figures 9 and 10. Figures 9 (a) and (b) show how the velocity changes in time across the width and depth of the channel respectively; the velocity oscillates and is zero at the walls as expected. For figure 9 (a), $z^{*}$ is chosen away from the walls to ensure non-zero velocity (and similarly for $y^{*}$ in figure 9 (b)). A time averaged velocity profile is plotted in figure 9 (c), this is qualitatively similar to the steady velocity profile in figure 2 (a) but the magnitude of the velocity is reduced.

The orientation parameter is shown over time and for multiple values of $z^{*}$ in figures 10 (a), for $y^{*}=1.9\times 10^{-3}~{}$m, and 10 (b), for $y^{*}=3.5\times 10^{-4}~{}$m. Towards the middle width (figure 10 (a)), the orientation parameter increases with the channel depth and in the centre of the channel (figure 10 (b)) the opposite occurs. Figure 10 (c) displays the time averaged orientation parameter which agrees well with that produced in steady flow (figure 2 (b)); this is also seen for the linear dichroism signal in figure 10 (d). Note here that the peak values of the orientation parameter are similar for steady and oscillatory flow.

Next, the spatially averaged orientation parameter $\overline{S}$ is compared for various pressure gradient and frequencies of oscillation (figure 11). As the pressure gradient $G^{*}$ increases, the orientation parameter increases; there is an indication that the impact of $G^{*}$ on $\overline{S}$ is reducing as $G^{*}$ becomes large. The impact of the frequency of oscillation $\omega^{*}$ is to decrease $\overline{S}$, a result seen throughout all figures. The relationship between the pressure gradient and frequencies of oscillation, for a choice of channel depth, is summarised in table 2. Changes in the number density have not been included for oscillating flow as the impact on the orientation parameter is small and reduces further for increasing frequencies of oscillations.

The orientation parameter is calculated for increasing channel depth $h^{*}$ and for a number of frequencies of oscillation and channel width (figure 12). Increasing $W^{*}$, shown by line markers in figure 12 (a), decreases the orientation parameter in general. The orientation parameter has a non-monotonic response to $h^{*}$as the frequency of oscillations is increased; for small $\omega^{*}$, increasing $h^{*}$ increases the orientation parameter whereas for large $\omega^{*}$ this relationship begins to reverse. This trend is highlighted in figure 12 (b), where $W^{*}=2\times 10^{-3}~{}$m. The relationship between $h^{*}$, $W^{*}$, $\omega^{*}$ and the orientation parameter $\overline{S}$ has been summarised in table 3; when $\omega^{*}$ is large, a larger reduction in $\overline{S}$ is seen with increasing channel width.

To further investigate these properties, the orientation parameter is plotted for increasing channel width for a range of channel depth and frequencies of oscillations (figure 13). As seen previously, the orientation parameter $\overline{S}$ increases with channel depth and decreases with increasing channel width. Increasing the frequency of oscillations reduces the orientation parameter in general; for larger channel depths, the decrease in the orientation parameter with $\omega^{*}$ is greater.

6 Discussion

The coupled pressure-driven flow and orientation dynamics of a dilute suspension of elongated particles has been modelled, to provide a mechanistic underpinning for flow linear dichroism spectroscopy. The model was applied specifically to a prototype hand-held device for waterborne pathogen detection, to explore the effect of changing the channel dimensions and pressure gradient, and to assess the feasibility of oscillatory flow. Mathematically, the problem involved coupling the Navier-Stokes equations, simplified via lubrication theory, to the Fokker-Planck equation for the shear-induced orientation dynamics. Whereas the flow problem involved two spatial independent variables and time, the orientation problem involved two angular independent variables and time, applied at each spatial location in the domain; the latter follows from the assumption that the particles are sufficiently small that they encounter a locally uniform shear rate. The flow and orientation problems were coupled via corrections to the fluid constitutive law involving moments of the orientation distribution, as described by established orientation Brownian suspension mechanics theory.

The system was solved numerically via spherical harmonics for the orientation distribution and a finite difference method for the fluid dynamics, with iterative coupling; the oscillatory problem was temporally discretised via a second-order explicit method for orientation dynamics and an alternating direction implicit method for the fluid dynamics. The oscillatory flow is rather more computationally expensive due to temporal dependence; a parameter set of 175 tuples was explored for the oscillatory flow, with each tuple requiring 14–20 hours of walltime on 5 workstation cores. The dependence of particle orientation on channel width, depth, pressure gradient, particle number density and frequency of oscillation were reported and compared for both the steady and oscillatory systems.

To model linear dichroism, the analysis focused on a single orientation parameter, calculated as the average of the difference between the components of the particle alignment parallel to, and perpendicular to, the flow direction. It was noted that the earlier (homogeneous shear) analysis of McLachlan [40] applied an alignment formula based on uni-axial orientation for elongated particles embedded in a membrane. However because flow linear dichroism produces a bi-axial orientation distribution – in other words there is no rotational symmetry about the flow direction in general – this earlier application was inappropriate. The correct bi-axial formula was compared with the uni-axial formula; it was found that the qualitative behaviour was similar, however the uni-axial definition overestimates the degree of orientation.

In both steady and oscillatory flows, increasing the pressure gradient increases the shear rate and hence the degree of alignment in the system. For steady flow, a plateau was observed close to a – perhaps maximum achievable – value of $\overline{S}=0.75$, above $G^{*}=5\times 10^{4}~{}$Pa/s and $h^{*}=2\times 10^{-3}~{}$m, whereas for oscillating flow the mean orientation did not exceed $0.43$ for the range of parameters analysed. This result does however suggest that oscillatory flow may be a viable method for producing practically-useful levels of alignment, which in turn may be valuable in the analysis of small sample volumes. Alignment decays monotonically with driving frequency, however an increase from $10$ rad/s to $40$ rad/s results in less than a 25% reduction in alignment. Increasing channel height has a significantly greater impact than increasing channel width, although for oscillating flow there is evidence of a local maximum in alignment as channel height is increased. For a given limited sample volume there is a clear trade-off between channel volume and oscillatory frequency; as channel volume is increased, a higher oscillatory frequency is needed to keep the sample underneath the detection window. The data reported should enable the optimal configuration for a given volume to be calculated.

Signal production was also modelled by multiplying the alignment parameter by the number density of particles; due to the relatively small impact of particle density on flow, the relationship between signal and number density was therefore approximately linear. While increasing the number density may therefore appear to be a simple way to optimise signal, in practice there are trade-offs with scattering – which was not modelled in the present case – and cost of the reagent. Optimising the remaining design parameters is therefore important.

Avenues for further work include experimental testing of the predictions of how signal should vary with channel dimensions and pressure gradient, and direct testing of oscillating flow linear dichroism spectroscopy. It will also be important to take account of scattering and hence produce a more accurate model of how number density affects signal production. Other avenues include incorporating flexibility of particles – since the persistence length of the M13 bacteriophage considered here is approximately $1.3$ $\mu$m [31] and the fibres are of length $0.8$ $\mu$m a rigid fibre approximation is reasonable, however future work could incorporate the flexibility of the fibres.

While the present manuscript focused on flows of passive particle suspensions in thin rectangular channels, future work may include investigating the behaviour of an active suspension, for example involving self-motile elongated cells such as spermatozoa in pressure-driven flow or complex behaviours such as trapping in shear flow [4]. Active suspensions display collective behaviour [28, 33, 48, 50] and superfluidity [36]. The model considered here can be augmented by an extra stress term to account for self-motile behaviour [44],

[TABLE]

The dynamics of suspensions of elongated particles provides an enduring subject in fluid dynamics, with novel application areas in continuing to be discovered. The challenging task of solving the complex models associated with non-homogeneous shear and time-dependent forcing is now becoming within reach of multicore computing hardware. We hope that the present manuscript provides both useful information in the design and optimisation of the technological system considered, and an effective framework which can be applied and extended to other systems in biophysical spectroscopy and cell analysis.

Funding

This work was supported by Biotechnology and Biological Sciences Research Council Industrial CASE Studentship (BB/L015587/1) and Engineering and Physical Sciences Research Council Healthcare Technologies Challenge Award (Rapid Sperm Capture EP/N021096/1).

Acknowledgements

The authors acknowledge Linear Diagnostics Limited for supporting this research. We thank Professor Tim Dafforn, Professor Alison Rodger, Dr Craig Holloway and Dr Meurig Gallagher for valuable discussion and help throughout this project.