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

G. Cupples  School of Mathematics, University of Birmingham, B15 2TT, U.K. D.
J. Smith  School of Mathematics, University of Birmingham, B15 2TT, U.K. M. R.
Hicks  Linear Diagnostics Ltd, BioHub Birmingham, 97 Vincent Drive,
Birmingham, B15 2SQ, U.K.  R. J. Dyson  School of Mathematics, University of
Birmingham, B15 2TT, U.K. Correspondence email<EMAIL_ADDRESS>
###### 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.

##  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.

Figure 1: Schematic of the flow of a fibre-laden fluid, forced by pressure
gradient $G^{*}$ with a distribution of orientation, in a thin rectangular
channel. (a) Rectangular channel of width $2W^{*}$ and depth $2h^{*}$ and the
$x^{*}$–axis is the molecular orientation axis. The orientation parameter is
calculated as the proportional difference in orientation parallel and
perpendicular to the molecular orientation axis, highlighted by the arrows to
the left of the channel. Note that for oscillating flow, the pressure gradient
is replaced with $G^{*}e^{i\omega^{*}t}$ for frequency of oscillations
$\omega^{*}$. (b) The orientation angles $\theta\in[0,\pi]$ and
$\phi\in[0,2\pi)$ away from the coordinate axis describe the orientation of
the fibres.

##  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

$\displaystyle\rho^{*}\left(\mathinner{\dfrac{\partial{}{\boldsymbol{u}^{*}}}{\partial{{t^{*}}}}}+(\boldsymbol{u}^{*}\cdot\nabla^{*})\boldsymbol{u}^{*}\right)$
$\displaystyle=\nabla^{*}\cdot\boldsymbol{\sigma}^{*},$ (1)
$\displaystyle\nabla^{*}\cdot\boldsymbol{u}^{*}$ $\displaystyle=0,$ (2)

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,

$\boldsymbol{u}^{*}(x^{*},\pm
W^{*},z^{*},t^{*})=0,\quad\boldsymbol{u}^{*}(x^{*},y^{*},\pm h^{*},t^{*})=0.$
(3)

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,

$\mathinner{\dfrac{\partial{}\psi}{\partial{{t^{*}}}}}+\nabla_{\boldsymbol{p}}\cdot(\boldsymbol{\Omega}^{*}\psi)=D_{r}^{*}\nabla_{\boldsymbol{p}}^{2}\psi,$
(4)

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

$\boldsymbol{p}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta),$ (5)

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

$\boldsymbol{\Omega}^{*}=(\textbf{I}-\boldsymbol{pp})\cdot(\alpha_{0}\boldsymbol{e}^{*}+\boldsymbol{\omega}^{*})\cdot\boldsymbol{p},$
(6)

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,

$\displaystyle\psi(\phi,\theta)$
$\displaystyle=\psi(\phi+\pi,\theta)=\psi(\phi,\theta+2\pi),$ (7)
$\displaystyle\int_{s}\psi d\boldsymbol{p}=1,$ (8)

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

$\boldsymbol{\sigma}^{*}=\boldsymbol{\sigma}_{I}^{*}+\boldsymbol{\sigma}_{D}^{*}+\boldsymbol{\sigma}_{P}^{*},$
(9)

where the Newtonian contribution is

$\boldsymbol{\sigma}_{I}^{*}=-P^{*}\textbf{I}+2\mu^{*}\boldsymbol{e}^{*},$
(10)

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

$\displaystyle\boldsymbol{\sigma}_{D}^{*}=$ $\displaystyle 2\mu^{*}\Phi
D_{r}^{*}\alpha_{r}\int_{s}\left(\boldsymbol{pp}-\frac{\textbf{I}}{3}\right)\psi\operatorname{d\\!}\boldsymbol{p},$
(11)
$\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}$
(12)

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

$\displaystyle\begin{split}u^{*}=\frac{G^{*}{h^{*}}^{2}}{\mu^{*}}u,\quad
v^{*}=\frac{\delta G^{*}{h^{*}}^{2}}{\mu^{*}}v,\quad w^{*}=\frac{\delta
G^{*}{h^{*}}^{2}}{\mu^{*}}w,\\\ \quad x^{*}=\frac{h^{*}}{\delta}x,\quad
y^{*}=h^{*}y,\quad z^{*}=h^{*}z,\quad
P^{*}=\frac{G^{*}h^{*}}{\delta}P,\end{split}$ (13)

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

$\displaystyle\begin{split}S=&\langle p_{x^{*}}^{2}\rangle-\langle
p_{y^{*}}^{2}\rangle,\\\
=&\langle\sin^{2}\theta\cos^{2}\phi\rangle-\langle\sin^{2}\theta\sin^{2}\phi\rangle,\end{split}$
(14)

for a 3D suspension where

$\langle\bullet\rangle=\int_{0}^{2\pi}\int_{0}^{\pi}\bullet\,\psi\sin\theta\operatorname{d\\!}\theta\operatorname{d\\!}\phi.$
(15)

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

$S=\frac{1}{2}\left(3\int_{s}\cos^{2}\theta_{S}\,\psi\operatorname{d\\!}\boldsymbol{p}-1\right),$
(16)

where

$\theta_{S}=\arccos(|\sin\theta\cos\phi|).$

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

$\frac{1}{6}\Lambda(\psi)=P_{e}\Upsilon(\psi),$ (17)

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

$\displaystyle\Lambda(\psi)$
$\displaystyle=\dfrac{1}{\sin{\theta^{\prime}}}\mathinner{\dfrac{\partial{}}{\partial{{\theta^{\prime}}}}}\left(\sin{\theta^{\prime}}\mathinner{\dfrac{\partial{}\psi}{\partial{{\theta^{\prime}}}}}\right)+\dfrac{1}{\sin^{2}{\theta^{\prime}}}\mathinner{\dfrac{\partial{{}^{2}}\psi}{\partial{{\phi^{\prime}}^{2}}}},$
(18) $\displaystyle\Upsilon(\psi)$
$\displaystyle=\dfrac{\sin{\phi^{\prime}}\cos{\phi^{\prime}}}{\sin{\theta^{\prime}}}\mathinner{\dfrac{\partial{}}{\partial{{\theta^{\prime}}}}}\left(\sin^{2}{\theta^{\prime}}\cos{\theta^{\prime}}\,\psi\right)-\mathinner{\dfrac{\partial{}}{\partial{{\phi^{\prime}}}}}\left(\sin^{2}{\phi^{\prime}}\,\psi\right),$
(19)

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

$\psi({\phi^{\prime}},{\theta^{\prime}})=\sum_{n=0}^{N}\sum_{m=0}^{n}\left(A_{0\,n}^{m}P_{n}^{m}(\cos{\theta^{\prime}})\cos(m{\phi^{\prime}})+A_{1\,n}^{m}P_{n}^{m}(\cos{\theta^{\prime}})\sin(m{\phi^{\prime}})\right),$
(20)

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

$\displaystyle\Lambda(P_{n}^{m}(\cos{\theta^{\prime}})\cos(m{\phi^{\prime}}))$
$\displaystyle=-n(n+1)P_{n}^{m}(\cos\theta^{\prime})\cos(m\phi^{\prime}),$
(21)
$\displaystyle\Lambda(P_{n}^{m}(\cos{\theta^{\prime}})\sin(m{\phi^{\prime}}))$
$\displaystyle=-n(n+1)P_{n}^{m}(\cos\theta^{\prime})\sin(m\phi^{\prime}),$
(22)
$\displaystyle\Upsilon(P_{n}^{m}(\cos{\theta^{\prime}})\cos(m{\phi^{\prime}}))$
$\displaystyle=-\sum_{p=m-2}^{m+2}\sum_{q=n-2}^{n+2}a_{n,q}^{m,p}P_{n}^{m}(\cos{\theta^{\prime}})\sin(m{\phi^{\prime}}),\quad
m\geq 0,$ (23)
$\displaystyle\Upsilon(P_{n}^{m}(\cos{\theta^{\prime}})\sin(m{\phi^{\prime}}))$
$\displaystyle=\sum_{p=m-2}^{m+2}\sum_{q=n-2}^{n+2}a_{n,q}^{m,p}P_{n}^{m}(\cos{\theta^{\prime}})\cos(m{\phi^{\prime}}),\quad
m>0.$ (24)

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

$\displaystyle\frac{q(q+1)}{6}A_{0\,q}^{p}$
$\displaystyle=-P_{e}\sum_{n=0}^{N}\sum_{m=0}^{n}a_{n,q}^{m,p}A_{1\,n}^{m},$
(25) $\displaystyle\frac{q(q+1)}{6}A_{1\,q}^{p}$
$\displaystyle=P_{e}\sum_{n=0}^{N}\sum_{m=0}^{n}a_{n,q}^{m,p}A_{0\,n}^{m},$
(26)

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

$\displaystyle\begin{split}-1=&\mathinner{\dfrac{\partial{}}{\partial{y}}}\left(\left\\{1+4\Phi\left[\alpha_{2}\int_{s}p_{1}^{2}p_{2}^{2}\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{3}}{2}\int_{s}(p_{1}^{2}+p_{2}^{2})\psi\,\operatorname{d\\!}\boldsymbol{p}+\frac{\alpha_{4}}{2}\right]\right\\}\mathinner{\dfrac{\partial{}u}{\partial{y}}}\right.\\\
&\hskip
102.43008pt\left.+4\Phi\left[\alpha_{2}\int_{s}p_{1}^{2}p_{2}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{3}}{2}\int_{s}p_{2}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}\right]\mathinner{\dfrac{\partial{}u}{\partial{z}}}\right)\\\
+&\mathinner{\dfrac{\partial{}}{\partial{z}}}\left(\left\\{1+4\Phi\left[\alpha_{2}\int_{s}p_{1}^{2}p_{3}^{2}\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{3}}{2}\int_{s}(p_{1}^{2}+p_{3}^{2})\psi\,\operatorname{d\\!}\boldsymbol{p}+\frac{\alpha_{4}}{2}\right]\right\\}\mathinner{\dfrac{\partial{}u}{\partial{z}}}\right.\\\
&\hskip
102.43008pt\left.+4\Phi\left[\alpha_{2}\int_{s}p_{1}^{2}p_{2}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{3}}{2}\int_{s}p_{2}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}\right]\mathinner{\dfrac{\partial{}u}{\partial{y}}}\right)\\\
+&\dfrac{2\Phi\alpha_{r}}{P_{G}}\left(\mathinner{\dfrac{\partial{}}{\partial{y}}}\int_{s}p_{1}p_{2}\psi\,\operatorname{d\\!}\boldsymbol{p}+\mathinner{\dfrac{\partial{}}{\partial{z}}}\int_{s}p_{1}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}\right),\end{split}$
(27)

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

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

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

$\mathinner{\dfrac{\partial{}{P^{*}}}{\partial{{x^{*}}}}}=G^{*}\exp(i\,\omega^{*}t^{*}),$
(29)

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

$\mathinner{\dfrac{\partial{}\psi}{\partial{\tau}}}=\frac{1}{6}\Lambda(\psi)-P_{e}(\tau)\Upsilon(\psi),$
(30)

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,

$\psi(\phi^{\prime},\theta^{\prime},t)=\sum_{n=0}^{N}\sum_{m=0}^{n}\left(A_{0\,n}^{m}(t)P_{n}^{m}(\cos\theta^{\prime})\cos(m\phi^{\prime})+A_{1\,n}^{m}(t)P_{n}^{m}(\cos\theta^{\prime})\sin(m\phi^{\prime})\right).$
(31)

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

$\displaystyle\mathinner{\dfrac{\operatorname{d\\!}{}A_{0\,q}^{p}}{\operatorname{d\\!}{t}}}$
$\displaystyle=-\dfrac{q(q+1)}{6}A_{0\,q}^{p}-P_{e}(t)\displaystyle\sum_{n=0}^{N}\sum_{m=0}^{n}a_{n,q}^{m,p}A_{1\,n}^{m},$
(32)
$\displaystyle\mathinner{\dfrac{\operatorname{d\\!}{}A_{1\,q}^{p}}{\operatorname{d\\!}{t}}}$
$\displaystyle=-\dfrac{q(q+1)}{6}A_{1\,q}^{p}+P_{e}(t)\displaystyle\sum_{n=0}^{N}\sum_{m=0}^{n}a_{n,q}^{m,p}A_{0\,n}^{m},$
(33)

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

$\displaystyle\begin{split}\alpha^{2}\mathinner{\dfrac{\partial{}u}{\partial{t}}}=&-\exp(it)+\dfrac{2\Phi\alpha_{r}}{P_{G}}\left(\mathinner{\dfrac{\partial{}}{\partial{y}}}\int_{s}p_{1}p_{2}\psi\,\operatorname{d\\!}\boldsymbol{p}+\mathinner{\dfrac{\partial{}}{\partial{z}}}\int_{s}p_{1}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}\right)\\\
+&\mathinner{\dfrac{\partial{}}{\partial{y}}}\left(\left\\{1+4\Phi\left[\alpha_{2}\int_{s}p_{1}^{2}p_{2}^{2}\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{3}}{2}\int_{s}(p_{1}^{2}+p_{2}^{2})\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{4}}{2}\right]\right\\}\mathinner{\dfrac{\partial{}u}{\partial{y}}}\right.\\\
&\hskip
71.13188pt\left.+4\Phi\left[\alpha_{2}\int_{s}p_{1}^{2}p_{2}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{3}}{2}\int_{s}p_{2}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}\right]\mathinner{\dfrac{\partial{}u}{\partial{z}}}\right)\\\
+&\mathinner{\dfrac{\partial{}}{\partial{z}}}\left(\left\\{1+4\Phi\left[\alpha_{2}\int_{s}p_{1}^{2}p_{3}^{2}\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{3}}{2}\int_{s}(p_{1}^{2}+p_{3}^{2})\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{4}}{2}\right]\right\\}\mathinner{\dfrac{\partial{}u}{\partial{z}}}\right.\\\
&\hskip
71.13188pt\left.+4\Phi\left[\alpha_{2}\int_{s}p_{1}^{2}p_{2}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}+\dfrac{\alpha_{3}}{2}\int_{s}p_{2}p_{3}\psi\,\operatorname{d\\!}\boldsymbol{p}\right]\mathinner{\dfrac{\partial{}u}{\partial{y}}}\right),\end{split}$
(34)

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)).

(a) Variable parameters

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

(b) Fixed parameters and physical constants

Parameter | Description | Value | Units
---|---|---|---
$a^{*}$ | Particle length | $8\times 10^{-7}$ | m
$b^{*}$ | Particle width | $6\times 10^{-9}$ | m
$V_{c}^{*}$ | Particle volume | $1.21\times 10^{-21}$ |  m3
$T^{*}$ | Temperature | $295$ | K
$\mu^{*}$ | Viscosity (water) | $9.5\times 10^{-4}$ | Pa s
$\rho^{*}$ | Fluid density (water) | $1\times 10^{3}$ |  kg/m3
$D_{r}^{*}$ | Rotational diffusion coefficient | $40.65$ |  s-1
$m_{w}^{*}$ | Particle molecular weight | $2\times 10^{7}$ | g/mol
$\epsilon^{*}$ | Extinction coefficient | $0.38$ |  m2/g
$K_{B}^{*}$ | Boltzmann constant | $1.38\times 10^{-23}$ |  kg m2/K s2
$N_{A}^{*}$ | Avogadro constant | $6\times 10^{23}$ | phage/mol

(c) Dimensionless groups

Parameter | Description | Definition | Value
---|---|---|---
$\Phi$ | Volume fraction | $V_{c}^{*}n_{d}^{*}$ | $1.11\times 10^{-5}$
$P_{G}$ | Péclet number (global) | $G^{*}h^{*}/D_{r}^{*}\mu^{*}$ | $34.17$
$\alpha^{2}$ | Womersley number | $\omega^{*}{h^{*}}^{2}\rho^{*}/\mu^{*}$ | $18.95$

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.

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

$\overline{\,\bullet\,}=\frac{1}{4W^{*}h^{*}}\int_{-W^{*}}^{W^{*}}\int_{-h^{*}}^{h^{*}}\bullet\operatorname{d\\!}y^{*}\operatorname{d\\!}z^{*},\quad\overline{\,\bullet\,}^{y^{*}}=\frac{1}{4W^{*}}\int_{-W^{*}}^{W^{*}}\bullet\operatorname{d\\!}y^{*}.$
(35)

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

$LD=\frac{\epsilon^{*}m_{w}^{*}n_{d}^{*}}{N_{A}^{*}}S_{rms},$ (36)

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: Velocity field, orientation parameter and predicted linear dichroism
signal for steady flow and the standard parameter set. (a) Velocity profile,
(b) orientation parameter and (c) linear dichroism signal.

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.

Figure 3: The orientation distribution function $\psi$, in orientation space
$\psi\in[0,\,2\pi]$, $\theta\in[0,\,\pi)$, at different points in the cross
sectional face of the flow using the standard parameter set. (a) The shear
rate $\dot{\gamma}^{*}$ in the cross sectional face of the channel. The dashed
circles, with corresponding labels, show the points in the flow where $\psi$
has been evaluated: (b) $y^{*}=-0.8\times 10^{-3}~{}$m, $z^{*}=-1\times
10^{-4}~{}$m, (c) $y^{*}=-0.4\times 10^{-3}~{}$m, $z^{*}=2\times 10^{-4}~{}$m,
(d) $y^{*}=0.95\times 10^{-3}~{}$m, $z^{*}=2.75\times 10^{-4}~{}$m and (e)
$y^{*}=0~{}$m, $z^{*}=0~{}$m.

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$.

Figure 4: The difference between the Newtonian and suspension flow fields in
the cross-sectional face of the channel. (a) Depicts the actual difference
$u_{N}^{*}-u^{*}$ and (b) depicts the difference relative to the magnitude of
the Newtonian velocity $(u_{N}^{*}-u^{*})/u_{N}^{*}$.

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}$.

Figure 5: A comparison between the width averaged orientation parameter
calculated in the current work (solid line) and the definition used by
McLachlan [40] (dashed line), for the standard parameter set.

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.

Figure 6: The orientation parameter and linear dichroism signal for steady
flow with increasing pressure gradient number density. (a) The orientation
parameter $\overline{S}$ with increasing pressure gradient, where
$n_{d}^{*}=1\times 10^{18}~{}$phage/m3. (b) The linear dichrosim signal
$\overline{LD}$ with increasing pressure gradient and three different particle
number densities: $n_{d}^{*}=1\times 10^{16}~{}$phage/m3 (dashed line),
$n_{d}^{*}=1\times 10^{17}~{}$phage/m3 (dot-dashed line) and
$n_{d}^{*}=1\times 10^{18}~{}$phage/m3 (dotted line); the arrow denotes
increasing number density. The other parameters are $h^{*}=3\times
10^{-4}~{}$m and $W^{*}=1\times 10^{-3}~{}$m.

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)).

Figure 7: The orientation parameter changing with channel dimensions for
steady flow. (a) $\overline{S}$ for a range of channel width and depth, (b)
$\overline{S}$ changing with increasing channel width $W^{*}$ for multiple
small channel width values: $h^{*}=2\times 10^{-4}~{}$m (solid line),
$h^{*}=3\times 10^{-4}~{}$m (dashed line), $h^{*}=4\times 10^{-4}~{}$m (dot-
dashed line) and $h^{*}=5\times 10^{-4}~{}$m (dotted line) and the arrow
denotes increasing channel width. Here $G^{*}=4.4\times 10^{3}~{}$Pa/m and
$n_{d}^{*}=7.33\times 10^{17}~{}$phage/m3.

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$.

Figure 8: Limit of $\overline{S}$ with increasing channel depth and pressure
gradient for steady flow. (a) $\overline{S}$ with varying channel depth
$h^{*}$ for a range of pressure gradients: $G^{*}=5\times 10^{3}~{}$Pa/m
(solid line), $G^{*}=1\times 10^{4}~{}$Pa/m (dashed line), $G^{*}=5\times
10^{4}~{}$Pa/m (dot-dashed line) and $G^{*}=10\times 10^{4}~{}$Pa/m (dotted
line). The arrow denotes increasing $G^{*}$. (b) $\overline{S}$ varying with
pressure gradient $G^{*}$ for a range of channel depths: $h^{*}=2\times
10^{-4}~{}$m (solid line), $h^{*}=6\times 10^{-4}~{}$m (dashed line),
$h^{*}=1\times 10^{-3}~{}$m (dot-dashed line) and $h^{*}=4\times 10^{-3}~{}m$
(dotted line). The arrow denotes increasing $h^{*}$.

###  5.2 Oscillating pressure gradient

Figure 9: Velocity field for oscillating flow and the standard parameter set.
(a) The velocity changing over time for multiple $y^{*}$ where $z^{*}=3\times
10^{-4}~{}$m. The arrows denote increasing positive values of $y^{*}$. (b) The
velocity changing over time for multiple points in $z^{*}$, where
$y^{*}=2\times 10^{-3}~{}$m. The arrows denote increasing positive values of
$z^{*}$. (c) The velocity field, averaged over half a time period, through the
cross-sectional face of the channel.

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.

Figure 10: The orientation parameter for oscillating flow and the standard parameter set. (a) $S$ evolving in time for different values of $z^{*}$ where $y^{*}=1.9\times 10^{-3}~{}$m. The arrow denotes increasing $z^{*}$. (b) $S$ evolving in time for different values of $z^{*}$ where $y^{*}=3.5\times 10^{-4}~{}$m. The arrow denotes increasing $z^{*}$. (c) $S$ through the cross-sectional face of the channel, averaged over half a time period. (d) $LD$ averaged over half a time period. $G^{*}$ | $h^{*}$ | $\overline{S}$ | $\overline{S}$
---|---|---|---
|  | ($\omega^{*}=10~{}$s-1) | ($\omega^{*}=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 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\times 10^{18}~{}$phage/m3 and
$W^{*}=2\times 10^{-3}~{}$m.

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.

Figure 11: Spatially-averaged orientation parameter, for oscillating flow,
varying with increasing pressure gradient and for three values of
$\omega^{*}$: $\omega^{*}=10~{}$s-1 (solid line), $\omega^{*}=20~{}$s-1
(dashed line) and $\omega^{*}=30~{}$s-1 (dotted line). Here $h^{*}=3\times
10^{-4}~{}$m, $W^{*}=1\times 10^{-3}~{}$m, $n_{d}^{*}=1\times
10^{17}~{}$phage/m3 and the arrow denotes increasing frequencies of
oscillation. Figure 12: Spatially averaged orientation parameter for
oscillating flow, changing with channel dimensions and frequency of
oscillation. (a) $\overline{S}$ varying with channel depth for multiple
frequencies of oscillation: $\omega^{*}=10~{}$s-1 (solid line),
$\omega^{*}=20~{}$s-1 (dashed line) and $\omega^{*}=30~{}$s-1 (dotted line).
The lines without markers are for $W^{*}=5\times 10^{-4}~{}$m and the lines
with markers are for $W^{*}=2\times 10^{-3}~{}$m. The arrow denotes increasing
frequencies of oscillations. (b) $\overline{S}$ changing with frequency of
oscillations for multiple values of $h^{*}$: $h^{*}=2\times 10^{-4}~{}$m
(solid line), $h^{*}=3\times 10^{-4}~{}$m (dashed line), $h^{*}=4\times
10^{-4}~{}$m (dot-dashed line) and $h^{*}=5\times 10^{-4}~{}$m (dotted line)
and for $W^{*}=2\times 10^{-3}~{}$m. For both figures, $G^{*}=8\times
10^{3}~{}$Pa/m and $n_{d}^{*}=1\times 10^{17}~{}$phage/m3.

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.

Figure 13: Spatially averaged orientation parameter, for oscillating flow,
changing with increasing channel width $W^{*}$. Three values of the frequency
of oscillations are plotted: $\omega^{*}=10~{}$s-1 (solid lines),
$\omega^{*}=20~{}$s-1 (dashed lines) and $\omega^{*}=30~{}$s-1 (dotted lines).
Two values of $h^{*}$ are also displayed: $h^{*}=2\times 10^{-4}~{}$m is shown
by plain lines and $h^{*}=5\times 10^{-4}~{}m$ is shown by the line markers.
The arrows denote increasing frequencies of oscillation.

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.

$W^{*}$ | $h^{*}$ | $\overline{S}$ | $\overline{S}$
---|---|---|---
|  | ($\omega^{*}=10~{}$s-1) | ($\omega^{*}=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$
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\times 10^{3}~{}$Pa/m and $n_{d}^{*}=1\times
10^{18}~{}$phage/m3.

##  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],

$\boldsymbol{\sigma}_{A}^{*}=\alpha_{1}\int_{s}\left(\boldsymbol{pp}-\frac{\textbf{I}}{3}\right)\psi\operatorname{d\\!}\boldsymbol{p}.$
(37)

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.

## References

  * [1] KL Aubrey and GJ Thomas Jr.  Raman spectroscopy of filamentous bacteriophage Ff (fd, M13, f1) incorporating specifically-deuterated alanine and tryptophan side chains. Assignments and structural interpretation.  Biophys J, 60(6):1337–1349, 1991.
  * [2] GK Batchelor.  The stress system in a suspension of force-free particles.  J Fluid Mech, 41(3):545–570, 1970.
  * [3] GK Batchelor.  The stress generated in a non-dilute suspension of elongated particles by pure straining motion.  J Fluid Mech, 46(4):813–829, 1971.
  * [4] RN Bearon and AL Hazel.  The trapping in high–shear regions of slender bacteria undergoing chemotaxis in a channel.  J Fluid Mech, 771, 2015.
  * [5] RB Bird, RC Armstrong, O Hassager, and CF Curtiss.  Dynamics of polymeric liquids, volume 1.  Wiley New York, 1977.
  * [6] BM Bulheller, A Rodger, and JD Hirst.  Circular and linear dichroism of proteins.  Phys Chem Chem Phys, 9(17):2020–2035, 2007.
  * [7] A Creppy, F Plouraboué, O Praud, X Druart, S Cazin, H Yu, and P Degond.  Symmetry-breaking phase transitions in highly concentrated semen.  J Royal Soc Interface, 13(123), 2016.
  * [8] G Cupples, RJ Dyson, and DJ Smith.  Viscous propulsion in active transversely isotropic media.  J Fluid Mech, 812:501–524, 2017.
  * [9] G Cupples, RJ Dyson, and DJ Smith.  On viscous propulsion in active transversely isotropic media.  J Fluid Mech, 855:408–420, 2018.
  * [10] T Daviter, N Chmel, and A Rodger.  Circular and linear dichroism spectroscopy for the study of protein–ligand interactions.  In Protein-Ligand Interactions, pages 211–241. Springer, 2013.
  * [11] M Doi and SF Edwards.  Dynamics of rod–like macromolecules in concentrated solution. Part 1\.  J Chem Soc Farad Trans 2, 74:560–570, 1978.
  * [12] J Douglas and DW Peaceman.  Numerical solution of two-dimensional heat-flow problems.  AIChE J, 1(4):505–512, 1955.
  * [13] RJ Dyson, JEF Green, JP Whiteley, and HM Byrne.  An investigation of the influence of extracellular matrix anisotropy and cell–matrix interactions on tissue architecture.  J Math Biol, 72(7):1775–1809, 2016.
  * [14] RJ Dyson and OE Jensen.  A fibre-reinforced fluid model of anisotropic plant cell growth.  J Fluid Mech, 655:472–503, 2010.
  * [15] JL Ericksen.  Transversely isotropic fluids.  Colloid Polym Sci, 173(2):117–122, 1960.
  * [16] B Ezhilan and D Saintillan.  Transport of a dilute active suspension in pressure–driven channel flow.  J Fluid Mech, 777:482–522, 2015.
  * [17] G Garab and H van Amerongen.  Linear dichroism and circular dichroism in photosynthesis research.  Photosyn Res, 101(2-3):135–146, 2009.
  * [18] MJ Glucksman, S Bhattacharjee, and L Makowski.  Three-dimensional structure of a cloning vector: X-ray diffraction studies of filamentous bacteriophage M13 at 7 å resolution.  J Mol Biol, 226(2):455–470, 1992.
  * [19] JEF Green and A Friedman.  The extensional flow of a thin sheet of incompressible, transversely isotropic fluid.  Eur J Appl Math, 19(03):225–257, 2008.
  * [20] E Guazzelli and JF Morris.  A physical introduction to suspension dynamics.  Cambridge University Press, 2011.
  * [21] NA Hill and TJ Pedley.  Bioconvection.  Fluid Dyn Res, 37(1):1–20, 2005.
  * [22] EJ Hinch and LG Leal.  Time–dependent shear flows of a suspension of particles with weak Brownian rotations.  J Fluid Mech, 57(4):753–767, 1973.
  * [23] EJ Hinch and LG Leal.  Constitutive equations in suspension mechanics. Part 1. General formulation.  J Fluid Mech, 71(03):481–495, 1975.
  * [24] EJ Hinch and LG Leal.  Constitutive equations in suspension mechanics. Part 2. Approximate forms for a suspension of rigid particles affected by Brownian rotations.  J Fluid Mech, 76(01):187–208, 1976.
  * [25] CR Holloway, G Cupples, DJ Smith, JEF Green, RJ Clarke, and RJ Dyson.  Influences of transversely isotropic rheology and translational diffusion on the stability of active suspensions.  Roy Soc Open Sci, 5(8):180456, 2018.
  * [26] CR Holloway, RJ Dyson, and DJ Smith.  Linear Taylor–Couette stability of a transversely isotropic fluid.  Proc R Soc Lond A, 471(2178), 2015.
  * [27] CR Holloway, DJ Smith, and RJ Dyson.  Linear Rayleigh–Benard stability of a transversely isotropic fluid.  Eur J Appl Math, pages 1–23, 2018.
  * [28] Y Hwang and TJ Pedley.  Bioconvection under uniform shear: linear stability analysis.  J Fluid Mech, 738:522–562, 2014.
  * [29] GB Jeffery.  The motion of ellipsoidal particles immersed in a viscous fluid.  Proc R Soc Lond A, 102(715):161–179, 1922.
  * [30] MR Kamal and AT Mutel.  The prediction of flow and orientation behavior of short fiber reinforced melts in simple flow systems.  Polym Composite, 10(5):337–343, 1989.
  * [31] A.S. Khalil, J.M. Ferrer, R.R. Brau, S.T. Kottmann, C.J. Noren, M.J. Lang, and A.M. Belcher.  Single m13 bacteriophage tethering and stretching.  Proc Natl Acad Sci, 104(12):4892–4897, 2007.
  * [32] S Kim and SJ Karrila.  Microhydrodynamics: principles and selected applications.  Courier Corporation, 2013.
  * [33] DL Koch and G Subramanian.  Collective hydrodynamics of swimming microorganisms: Living fluids.  Annu Rev Fluid Mech, 43:637–659, 2011.
  * [34] LG Leal and EJ Hinch.  The rheology of a suspension of nearly spherical particles subject to Brownian rotations.  J Fluid Mech, 55(04):745–765, 1972.
  * [35] FM Leslie.  Some constitutive equations for liquid crystals.  Arc Rational Mech Anal, 28(4):265–283, 1968.
  * [36] HM López, J Gachelin, C Douarche, H Auradou, and E Clément.  Turning bacteria suspensions into superfluids.  Phys Rev Lett, 115(2):028301, 2015.
  * [37] M Manhart.  Rheology of suspensions of rigid–rod like particles in turbulent channel flow.  J Non-Newton Fluid, 112(2):269–293, 2003.
  * [38] C Marchioli, M Fantoni, and A Soldati.  Orientation, distribution, and deposition of elongated, inertial fibers in turbulent channel flow.  Phys Fluids, 22(3):033301, 2010.
  * [39] R Marrington, TR Dafforn, DJ Halsall, JI MacDonald, M Hicks, and A Rodger.  Validation of new microvolume Couette flow linear dichroism cells.  Analyst, 130(12):1608–1616, 2005.
  * [40] JRA McLachlan, DJ Smith, NP Chmel, and A Rodger.  Calculations of flow-induced orientation distributions for analysis of linear dichroism spectroscopy.  Soft Matter, 9(20):4977–4984, 2013.
  * [41] PH Mortensen, HI Andersson, JJJ Gillissen, and BJ Boersma.  Dynamics of prolate ellipsoidal particles in a turbulent channel flow.  Phys Fluids, 20(9):093302, 2008.
  * [42] Z Niu, MA Bruckman, B Harp, CM Mello, and Q Wang.  Bacteriophage M13 as a scaffold for preparing conductive polymeric composite fibers.  Nano Res, 1(3):235–241, 2008.
  * [43] R Pacheco-Gómez, J Kraemer, S Stokoe, HJ England, CW Penn, E Stanley, A Rodger, J Ward, MR Hicks, and TR Dafforn.  Detection of pathogenic bacteria using a homogeneous immunoassay based on shear alignment of virus particles and linear dichroism.  Anal Chem, 84(1):91–97, 2011.
  * [44] TJ Pedley and JO Kessler.  A new continuum model for suspensions of gyrotactic micro-organisms.  J Fluid Mech, 212:155–182, 1990.
  * [45] A. Peterlin and H.A. Stuart.  Zur theorie der strömungsdoppelbrechung von kolloiden und großen molekülen in lösung.  Z Phys, 112(1-2):1–19, 1939.
  * [46] C Pozrikidis.  Introduction to theoretical and computational fluid dynamics.  Oxford University Press, 2011.
  * [47] A Rodger and B Nordén.  Circular dichroism and linear dichroism.  Oxford University Press, USA, 1997.
  * [48] D Saintillan and MJ Shelley.  Instabilities, pattern formation, and mixing in active suspensions.  Phys Fluids, 20(12):123304, 2008.
  * [49] SR Strand, S Kim, and SJ Karrila.  Computation of rheological properties of suspensions of rigid rods: stress growth after inception of steady shear flow.  J Non-Newton Fluid, 24(3):311–329, 1987.
  * [50] G Subramanian, DL Koch, and SR Fitzgibbon.  The stability of a homogeneous suspension of chemotactic bacteria.  Phys Fluids, 23(4):041901, 2011.


# A structural determinant of the behavior of water at hydration and
nanoconfinement conditions

Nicolás A. Loubet1 Alejandro R. Verde1 Gustavo A. Appignanesi1∗ 1 INQUISUR,
Departamento de Química, Universidad Nacional del Sur (UNS)-CONICET, Avenida
Alem 1253, 8000 Bahía Blanca, Argentina
* Corresponding author<EMAIL_ADDRESS>

###### Abstract

Liquid water is reluctant to lose hydrogen-bond coordination. Here we reveal
that it also demands contraction and reorientation of the second molecular
shell to compensate for coordination defects. Such molecular principle will be
shown to lie at the heart of the two-liquids scenario and to determine the
behavior of water in hydration and nanoconfinement. Furthermore, it allows to
define conditions for wettability (quantifying hydrophobicity and predicting
drying transitions), thus opening the possibility to elucidate the active role
of water in central fields of research.

Liquid water is not only a fascinating and challenging system in itself, but
its behavior is also determinant for a huge number of phenomena that are
central to our life and that impact in fields that range from biology to
materials science [1, 2, 3, 4, 5, 6, 7]. In turn, water is known to present
several structural, thermodynamic and dynamic anomalies which seem to stem
from its ability to coexist in two competing different liquid forms (low
density liquid, LDL, and high density liquid, HDL), as described by the two-
liquids scenario[1, 2, 3, 4, 5]. This description emerged from the proposition
of the existence of a second liquid-liquid critical point, which has recently
received strong validation[8, 4, 5, 9, 10, 11, 12, 13, 14, 15, 16, 17]. Some
of these anomalies are pivotal for making water irreplaceable by normal
liquids in contexts where its plays an active role. Remarkably, water is
considered to be an indispensable condition for the existence of life[1, 2];
indeed, when we search for the possibility of life outside planet earth, what
we do in fact is to look for water[1]. In turn, water’s relatively strong
intermolecular interactions, particularly its tetrahedrally-directed hydrogen
bonds, are essential in this anomalous behavior that distinguishes it from
simple liquids[1, 2, 3, 4, 5]. In hydration conditions, that is, when water is
in contact with a solute, a biomolecule or any other material, a fraction of
the water-water interactions should be replaced by the surface in order to
produce wetting[2, 18, 19, 20]. In some cases where this compensation is
provided by means of water-surface HBs or other significant water-surface
interactions, a good wetting is achieved, and the surface is termed as
hydrophilic[2, 21, 22]. However, the high reluctance of water to lose HB-
coordination[18, 19, 20] poses strong constraints to surfaces that interact
weakly with water so that this compensation cannot be fulfilled, thus
resulting in a poor wetting or hydrophobic behavior [2, 21, 22]. Nevertheless,
a full comprehension of the relative hydrophobicity of the different surfaces
in molecular terms is still lacking. An example in this regard is the
unexpected high hydrophilicity of graphene-like surfaces[23, 24, 25], which we
shall hereby consider. Additionally, the local geometry (as in the case of
cavities, pores, etc.) might conspire with the possibility for water to retain
a proper HB environment, even to the point to promote drying in certain
conditions[25]. Indeed, nanoconfinement affects dramatically some properties
of water and, if a proper HB network cannot be established, the liquid becomes
thermodynamically unstable[6] as, for example, in the case of the hydrophobic
collapse that promotes protein folding or in the case of biomolecular binding,
where the removal of labile water by a ligand has been established as a
principal driving force[18, 19, 20, 26, 27]. However, water indeed penetrates
certain subnanometric cavities, as it happens in small-radius nanotubes, a
fact long ago regarded as “surprising” on energetic grounds[28] and which we
shall hereby revisit. The above-expounded scenario reveals the fact that the
understanding of the principles of water hydration, including the development
of a theory of wettability and hydrophobicity at the nanoscale, is imperative
both from a fundamental and an applied standpoint. Notwithstanding, this
picture is still far from being complete [29, 6]. To move in such direction,
let us pose the following ansatz on the molecular structuring of the HB
network of liquid water: Provided its strong reluctance to lose its
tetrahedrally-directed coordination[18, 19, 20], anytime a water molecule
resigns a HB, certain directional energetic compensation would be demanded.
Such a tendency to retain intermolecular tetrahedrality results from the need
to avoid the otherwise uncompensated partial molecular charges. For bulk
water, we shall show that this requisite is met by a fine-tuned interplay
between the first and second neighbor shells, in line with the two-liquids
scenario. In turn, in the case of surface hydration and water nanoscale
confinement, this behavior should also be locally mimicked by the interacting
material in order to ensure hydrophilic wettability, as we shall also
demonstrate. Thus, this operational principle will be first validated in bulk
water by a simple structural metric which will then be applied to model
surfaces and nanoscale confinement to exhibit its potential to characterize
hydrophobicity and predict filling and drying transitions.

To begin, we introduce a new simple structural index to characterize water
structure in terms of directional (tetrahedral) intermolecular interactions.
We shall call this parameter $V_{4S}$, since it computes the different
energetic contributions (of other water molecules and also including the heavy
atoms of interacting surfaces for non-bulk situations) at four sites
tetrahedrally arranged around each water molecule (see Supplementary
Material). We sort the resulting energies in decreasing order so that $V_{4S}$
is the most weakly interacting site. In this regard, it is relevant to note
that previous structural indices for water are only appropriate for bulk water
and not suitable for other contexts where they would produce artifacts, as we
have indicated[30].

Figure 1: Distributions of the $V_{4S}$ index calculated at the inherent
dynamics for TIP3P water at 1 bar for T=300K. We also discriminate the
contributions of the T molecules (low local density, with good tetrahedral HB-
coordination) and D molecules (3-fold HB-coordinated and high local density).
Inset: an example of the environment of a D central molecule, where we
indicate in gray the molecules that form HB with it, in green the other
attractive ones (more than $k_{B}T$, where $k_{B}$ is Boltzmann’s constant),
in red the repulsive ones and in white the ones whose interaction with the
central one is almost null.

Fig. 1 displays the distribution of the $V_{4S}$ index, showing clear signs of
bimodality. Such behavior is even more evident if we additionally classify the
water molecules in terms of the previously defined $V_{4}$ index which, for a
given central molecule, finds the four water molecules with which it interacts
more strongly and picks the value of the less intense of such interactions[31,
32]. This parameter has been successful in dissecting water molecules into
LDL-like (called T molecules) and defects (D molecules)[31, 32]. T molecules
are characterized by a good four-fold tetrahedral HB-coordination in their
first shell (interactions stronger than 12 kJ/mol, compatible with good
quality HBs), while the D ones present three-fold HB coordination and are
responsible for the existence of the HDL state[31, 32]. Specifically, when we
discriminate the relative contributions of the T and D molecules to the
$V_{4S}$ distribution of Fig. 1, we find that the former ones peak to the left
(at around -17 kJ/mol), while the latter are broadly distributed to the right
with a mean value of $V_{4S}=-6.08$ kJ/mol, less intense but still quite
attractive. From one side, it is evident that T molecules leverage water’s
unusual ability to expand its structure to gain energy: an expansion of the
second neighbors shell allows the first one to enjoy a full 4-fold tetrahedral
HB-coordination with the central molecule. Thus, the $V_{4S}$ value for the T
molecules is dominated by the energy of a first-shell HB with virtually no
contribution of the second shell. D molecules, in turn, imply a broken HB,
since-second shell molecules have penetrated the intershell gap (thus
increasing the local density) and, hence, their first shell comprises only 3
HB-coordinated neighbors. However, as already recognized[33], the largely
overlooked HDL state is far from being unstructured. From the results of the
$V_{4S}$ distribution, it is now evident that the closer second-shell
molecules tend to arrange themselves in dipolar orientations to compensate the
partial molecular charge of the central molecule at the tetrahedral site of
the lacking HB (this can be achieved by many different local arrangements; see
the example shown in the inset of Fig. 1). In this sense, the central molecule
recovers part of the HB-loss, as evidenced by the peak at an attractive value
of -6.08 kJ/mol. We note that two or three second-shell molecules are usually
engaged in such an energetic compensation. Hence, while the expanded LDL state
is dominated by first-shell energetics, the high local density molecules
display an enhanced contribution of the partially contracted second-shell.
This interplay between structure and energetics (stemming form the directional
nature of the HB interactions) evidences the unique ability of water to
present, unlike simple liquids, a marked structural plasticity by taking
advantage of two opposed structural/energetic trends. Hence, a first
interesting lesson that emerges from these results is that the existence of
the LDL and HDL molecular states in bulk water might result from this fine-
tuning in the arrangement of the two first molecular shells in order to
satisfy the energetic constraints of the previously enunciated ansatz.

Figure 2: Distributions of the $V_{4S}$ index for TIP3P water at 1 bar for
T=300K close to the $SAM-OH$, the $SAM-CH_{3}$ and inside a small radius CNT.
We also include the case of bulk water.

In order to try to understand the structural basis of hydration, we decided to
investigate a couple of model systems: self-assembled monolayers (SAMs)
functionalized to act as hydrophobic or hydrophilic surfaces[22, 24, 25] and a
carbon nanotube, CNT (see Supplementary Material). In Fig. 2 we show results
of the $V_{4S}$ parameter calculated for the water molecules located closer
than 3.5Å from the SAMs in order to include the first peak of the water
density profile normal to the surface. In the case of the hydrophilic SAM
($SAM-OH$), the resulting distribution is almost identical to that of bulk
water (which we also include for comparison). In this case, the water
molecules tend to form three water-water HBs and an energetically analogous
fourth one is formed with the surface hydroxyls, thus creating an environment
equivalent to that of the 4-fold coordinated bulk molecules. For the
hydrophobic ($SAM-CH_{3}$) case, however, a water molecule close to the
surface is unable to retain full HB-coordination and, thus, is forced to
resign one HB in the direction of the surface. At such tetrahedral site the
$V_{4S}$ parameter shows that it feels a virtually null compensation from the
surface (thus producing a sharp peak close to zero). The crossing between this
peak and the distribution for the hydrophilic SAM (and also with the
distribution of $V_{4S}$ for bulk water) occurs at a value close to $k_{B}T$.
In fact, there is barely no population below $k_{B}T$ for the hydrophobic SAM.
These results reveal the ability of our new index to accurately discriminate
between hydrophobic and hydrophilic behaviors. For comparison, we note that an
appropriate measure of hydrophobicity is provided by the extent of water
density fluctuations and, thus, can be quantified by the probability of cavity
creation in small volumes, P(N=0), where N is the number of water molecules
inside the observation volume (usually spheres of radius r=3.3Å, compatible
with the size of a methane molecule, tangent to the surface)[22]. A $SAM-OH$
like the one we study here has been shown to present a value lower than
$P(N=0)=2\times 10^{-4}$ (low vacating probability and, thus, typical
hydrophilic behavior)[24]. This value is also almost identical to that in bulk
water[22]. In turn, the value for the $SAM-CH_{3}$ amounts to around $3\times
10^{-2}$ (a high dehydration propensity or, equivalently, a low work for water
removal or cavity creation, thus implying hydrophobicity)[24]. We also plot in
Fig. 2 the $V_{4S}$ distribution for the water molecules in contact with the
interior region of a carbon nanotube (we note that we also obtain a similar
result for a planar graphene sheet). The distribution of $V_{4S}$ for the CNT
displays a peak with a mean value of -5.99 kJ/mol well within the hydrophilic
region (although less negative than that for the $SAM-OH$). The hydrophilic
nature of graphene-based surfaces has been recognized both theoretically[24,
25] and by experimental water contact angle measurements over clean graphene-
based surfaces (that showed that the rather hydrophobic-like values reported
in previous studies had been in fact a result of aging and that arose from
hydrocarbon adsorption from the air)[23]. At this point we note that the
calculation of water vacating probabilities within surface-tangent small
radius spheres (P(N=0)) gave virtually identical values for graphene, CNTs of
different radius and $SAM-OH$ surfaces[24, 25]. This fact was quite
intriguing, since it implied that these surfaces presented an almost identical
dehydration propensity despite their different attraction intensities,
suggesting that the work for water removal would be the same once overcoming
certain threshold value and, thus, pointing to a saturation value for
hydrophilicity. Indeed, graphene-like surfaces present a high carbon density
and, thus, a neighboring water molecule would simultaneously interact with
many carbon atoms, adding up several van der Waals interactions. This fact
enables the establishment of full hydrophilic wetting, making evident that
nonpolar and hydrophobic are not synonyms. However, the total attraction that
a water molecule experiences with a graphene-like surface is still far less
intense as compared with a HB (roughly 35-40% of the value that water attains
at a $SAM-OH$ or in bulk conditions). Such a conundrum can be solved if we
refer to the behavior of the $V_{4S}$ index, since the mean value of this
indicator for the graphene-like surfaces is almost identical to the mean value
of the (high local density) D water molecules of the bulk water distribution
we have shown in Fig. 1 (we also recall that the distribution of $V_{4S}$ for
bulk water is, in turn, almost identical to the one for the $SAM-OH$). This
means that the water molecules in contact with a graphene-like surface receive
at their HB-uncoordinated tetrahedral site a compensation that amounts to the
value of the net interaction that a D molecule feels on average at its HB-
lacking site in bulk conditions($V_{4S}$ around -6kJ/mol). Thus, to attain
wettability and to reach hydrophilic-like behavior, it would not be necessary
to recover full HB stabilization, as it happens at the $SAM-OH$ (LDL-like
environments), but it would suffice to reach the local energetic requirement
of a typical D (HDL-like) molecule. In turn, if the interaction with any given
surface implies a $V_{4S}$ value larger than $V_{4S}$=-6kJ/mol, the surface
would exhibit a reduced wetting capability (reduced hydrophilicity), until
attaining a neat hydrophobic behavior when the $V_{4S}$ value reaches that of
the $SAM-CH_{3}$, close to -3kJ/mol (and wetting is unfavourable as compared
with a bulk D-state). Our results, thus, not only rationalize the
hydrophobicity scale for these different surfaces, but also extends the
validity of the two-liquids scenario to the hydration context by revealing the
conditions for wettability.

Figure 3: Top left: We show, for varying values of $\epsilon_{C-O}$, the mean
number of molecules inside the CNT (we exclude from the calculation the
regions at both CNT ends -2Å length each- in order to avoid for possible
border effects); top right: we depict the corresponding average values of
$V_{4S}$; bottom right: calculated distributions of $V_{4S}$ (solid lines) and
$V_{3S}$ (dashed lines) for chosen values of $\epsilon_{C-O}$; bottom left: we
include a detail of Fig. 2 for comparison. The red dashed lines indicate the
estimated values for complete filling and complete drying.

Concerning nanoconfinement, persistent penetration of small-diameter
(subnanometric) carbon nanotubes by a one-dimensional arrangement of water
molecules in simulations has long ago been deemed “surprising”[28], provided
the fact that only a fraction of the energetic loss (resignation of HB
coordination) could be recovered by water-wall attractions. Additionally, a
partial reduction of such attractions promoted drying[28]. To demonstrate that
such drying transitions can not only be explained by our approach but that
also their range of occurrence can be quantitatively predicted, we studied a
narrow CNT (see Supplementary material), varying the Lennard-Jones parameters
from $\epsilon_{C-O}=0.478369\ kJ~{}{mol}^{-1}$ and $\sigma_{C-O}=3.27510$ Å
to $\epsilon_{C-O}=0.228720\ kJ~{}{mol}^{-1}$ and $\sigma_{C-O}=3.44154$ Å,
including the situations studied in[28]. In Fig. 3 we provide the equilibrium
average number of water molecules inside the CNT as we vary the Lennard-Jones
parameters. Two evident regimes are detected close to each extreme $\epsilon$
values, as expected: completely filled for large $\epsilon$ and completely
desorbed at low $\epsilon$. These two regimes are separated by a transition
regime with a partial occupancy that results from an alternance of empty and
(totally or partially) filled periods. In turn, starting from many CNT-filled
configurations at $\epsilon_{C-O}=0.478369\ kJ~{}{mol}^{-1}$, we modified the
Lennard-Jones parameters within the indicated linear regime, and we calculated
the resulting $V_{4S}$ distributions. We also calculated the $V_{3}S$
distributions, since the water molecules inside the CNT resign two HBs (we
note, however, that the values of $V_{4S}$ and $V_{3}S$ are very similar, as
expected). By so doing, we generated plots relating $\epsilon_{C-O}$ to mean
$V_{4S}$ values, and we also included them in Fig. 3. When we compare these
plots, we find that starting from the highest $\epsilon_{C-O}$ value, the
completely filled regime of the CNT ends at a value of $\epsilon_{C-O}$ that
produces a mean $V_{4S}$ value very close to -6kJ/mol, the one characteristic
of the three-fold coordinated D water molecules in bulk conditions (cf. Fig.
1). In turn, when we reach an $\epsilon_{C-O}$ value that yields a $V_{4S}$
close to -3kJ/mol (around the crossing between the $SAM-OH$ and $SAM-CH_{3}$,
that is, the value which we have identified above as the beginning of the
fully hydrophobic regime), the CNT becomes to be completely desorbed. Thus, it
is evident that $V_{4S}$ is successful in quantitatively predicting
filling/drying transitions in this nanoconfinement setting. Again, complete
filling is compatible with a situation in which the CNT wall provides the
water molecules with an environment equivalent to that of the high local
density D molecules. It is not surprising, thus, that a full balance of the HB
loss was not required for water in order to penetrate the CNT, but that the
establishment of conditions compatible with the two-state scenario for water
is demanded. Below such value of $\epsilon_{C-O}$, filling begins to be only
partial, and the onset of complete drying occurs when the CNT walls are
equivalent to the $SAM-CH_{3}$, as the value of $V_{4S}$ falls close to
$k_{B}T$.

We have shown that while liquid water presents a strong tendency to maximize
the number of intermolecular hydrogen bonds (gaining energy by improving
4-fold coordination in the first-shell as the second neighbors are expanded),
HB-coordination defects demand certain local compensation of the molecular
charge at their HB-lacking position by second-shell contraction and
orientation. This structural-energetic interplay gives rise to two preferred
local molecular structures, as revealed by a new metric we have hereby
introduced, providing a detailed molecular rationalization of the two-liquids
scenario. Moreover, we also show that the same molecular principle is
operational at hydration and nanoconfinement conditions. In this sense, we
explain both the unexpected high hydrophilicity of graphene-like systems
(comparable to systems capable of stronger interactions with water and, thus,
pointing to a saturation hydrophilicity threshold) and the “surprising”
penetration of water in small radius carbon nanotubes (where the water-wall
interactions cannot fully recover the consequent loss in hydrogen bonding). In
both cases, it suffices that the surface compensates the tetrahedral HB-
lacking site at the level required by a high local density water molecule,
thus revealing a main molecular rationale underlying hydration, wetability and
nanoconfinement behavior for water and extending the validity of the two-
liquids scenario to these realms. Furthermore, the new metric we introduce to
test the molecular principle hereby revealed (which, at variance from the
situation of previous structural indicators is suitable for interfaces and
nanoconfinement) enables us to quantitatively characterize hydrophobicity and
to accurately predict water filling and drying transitions and, thus, could be
of great interest for different systems that operate in aqueous environments
within widespread contexts in biophysics and materials science.

##  I Supplementary Material: Methods

In this work we introduce a new simple structural index to characterize water
structure in terms of directional (tetrahedral) intermolecular interactions.
We shall call this parameter $V_{4S}$, since it computes the different
energetic contributions at four sites tetrahedrally arranged around each water
molecule. $V_{4S}$ is calculated for any given central water molecule as
follows: First, the H-O-H angle is opened until it reaches the value of a
tetrahedral angle, $arccos(-1/3)$ (this step is not necessary for models in
which the H-O-H angle is tetrahedral, like SPC/E). Subsequently, we define two
tetrahedral sites located at 1Å away from the oxygen in both directions and
geometrically complete the two other tetrahedrally oriented sites (also
located 1Å away from the oxygen). In this way, the four points/sites of a
perfect tetrahedron are determined. Thus, this method partitions space in four
regions as dictated by the tetrahedral arrangement of the four sites: We
select all the neighboring water molecules that are within a distance of 5Å
from each of the four sites. Then, we assign each of these molecules to their
closest site. We chose 5Å in order to include roughly both the first and
second neighbor shells of the central molecule in the case of bulk water. In
turn, for the situations in which we consider the interaction with a surface,
we also include all the heavy atoms of the surface that lie within the spheres
of radius 5Å centered at the four tetrahedral sites of the water molecule of
interest, and assign each of them to their closest site. Then, for each
tetrahedral site we add up the contributions of all the pair-wise interaction
potentials between the central molecule and all its assigned water molecules,
and also all its assigned heavy atoms in the case of surfaces, taking into
account Lennard-Jones and the Coulomb interactions. With this procedure, we
are left with four potential values, one for each point/site of the
tetrahedron. These are ordered so that $V_{1S}$ is the most interacting, and
$V_{4S}$ is the least interacting. As for other previous water structural
indicators like $V_{4}$[31, 32], both for bulk water and for the water
molecules in contact with the different systems studied, we calculated the
index for minimized configurations (inherent configurations) by using the
steepest-descent method to avoid the blurring effect of thermal energy[31,
32]. We note that this index is suitable for generic contexts, since the
energy contributions we consider at each of these four sites could arise
either from the interaction of the oxygen of the central molecule with its set
of neighboring water molecules and/or with the heavy atoms of any other
system. Thus, it could be applied both in bulk conditions or at interfaces or
nanoconfinement.

In order to try to understand the structural basis of hydration, we decided to
apply our new structural indicator to model systems: self-assembled monolayers
(SAMs) functionalized to act as hydrophobic or hydrophilic surfaces[22, 24,
25] and carbon nanotubes, CNT. Again, we note that previous structural
indicators for water would not be appropriate in this context[30]. Regarding
the simulations of the CNT, a 112-carbon (7.3 Å diameter) nanotube was
solvated in a box of 4554 TIP3P water molecules. Using the GROMACS package[34]
and the General AMBER force field (GAFF), this system was stabilized in NVT
conditions at 300K for several structural relaxation times ($\tau_{\alpha}$),
and subsequently stabilized in NPT conditions at 1bar and the same
temperature. Finally, a 40ns MD simulation was carried out, printing
configurations every 0.4ps. In all the stabilizations and dynamics, the
initial position of the nanotube was restrained (center of the box). The
interaction parameters used were $\epsilon_{C-O}=0.478369kJ~{}{mol}^{-1}$ and
$\sigma_{C-O}=3.27521$Å. In addition to this dynamics, 16 other 30ns dynamics
were performed each with configurations every 6ps, linearly varying the
interaction parameters from $\epsilon_{C-O}=0.478369kJ~{}{mol}^{-1}$ and
$\sigma_{C-O}=3.27510$Åto $\epsilon_{C-O}=0.228720kJ~{}{mol}^{-1}$ and
$\sigma_{C-O}=3.44154$Å. In the case of the SAMs[24, 25], we aligned 81 chains
($9\times 9$, chain separations of $4.35$ Å) of n-decane
($CH_{3}-(CH_{2})_{8}-CH_{3}$) in a parallel way so that their heads adopted a
square-like arrangement facing water (mimicking the monolayer structure
adopted by stearic acid at a water interface). The chain heads were
functionalized in two different ways: they were ended with methyl groups (for
the $SAM-CH_{3}$) or hydroxyls ($SAM-OH$). The system was solvated with 10654
TIP3P water molecules in an orthogonal cubic box with periodic boundary
conditions (that extended more than $20$ Å away from all the monolayer faces)
and was equilibrated within a canonical NVE ensemble with a Berendsen
thermostat [35]. All simulations were done using the AMBER10 molecular
simulation suite[36] with a 2 fs time step (we used the GAFF and FF99SB force
fields) and all calculations were performed in the NPT ensemble with a
Langevin thermostat.

## References

  * [1] P. Ball. Proc. Natl. Acad. Sci. U. S. A. 116, 13327 (2017).
  * [2] M. Chaplin. Nat Rev Mol Cell Biol 7, 861 (2006).
  * [3] P. Gallo, et al. Chem. Rev. 116, 7463 (2016).
  * [4] P. Gallo et al. Eur. Phys. J. E 44, 143 (2021).
  * [5] H. R. Corti, et al. Eur. Phys. J. 44 136 (2021).
  * [6] N. Giovambattista, P. J. Rossky and P. G. Debenedetti. Annu. Rev. Phys. Chem. 63, 179 (2012).
  * [7] R. Shi and H. Tanaka. Proc. Natl. Acad. Sci. U. S. A. 117, 26591 (2020).
  * [8] P. H. Poole, F. Sciortino, U. Essmann and H. E. Stanley. Nature 360, 324 (1992).
  * [9] F. Sciortino, I. Saika-Voivod and P. H. Poole. Phys. Chem. Chem. Phys. 13, 19759 (2011).
  * [10] J. C. Palmer, F. Martelli, Y. Liu, R. Car, A. Z. Panagiotopoulos and P. G. Debenedetti. Nature 510, 385 (2014).
  * [11] J. C. Palmer, P. H. Poole, F. Sciortino and P. G. Debenedetti. Chem. Rev. 118, 9129–9151 (2018).
  * [12] A. Nilsson and L. G. M. Pettersson. Chem. Phys. 389, 1 (2011).
  * [13] K. H. Kim, et al. Science 358, 1589 (2017).
  * [14] N. Stern, M. Seidl-Nigsch and T. Loerting. Proc. Natl. Acad. Sci. U. S. A. 116, 9191 (2019).
  * [15] K. H. Kim, et al.. Science 370, 978 (2020).
  * [16] J. W. Biddle, et al. J. Chem. Phys. 146, 034502 (2017).
  * [17] P. G. Debenedetti, F. Sciortino and G. H. Zerze. Science 369, 289 (2020).
  * [18] E. Schulz, M. Frechero, G. A. Appignanesi and A. Fernández. PLoS ONE 5, e12844 (2010).
  * [19] A. Fernández. Phys. Rev. Lett. 108, 188102 (2012).
  * [20] A. Fernández, Transformative Concepts for Drug Design: Target Wrapping. Springer-Verlag. (2010), ISBN-13 978-3642424588.
  * [21] N. Giovambattista, P. G. Debenedetti and P. J. Rossky. Proc. Natl. Acad. Sci. U. S. A. 106, 15181 (2009).
  * [22] S. N. Jamadagni, R. Godawat and S. Garde. Annu. Rev. Chem. Biomol. Eng. 2, 147 (2011).
  * [23] Z. Li, et al. Nature Materials 12, 925 (2013).
  * [24] S. R. Accordino, J. M. Montes de Oca, J. A. Rodriguez Fris and G. A. Appignanesi. J. Chem. Phys. 143, 154704 (2015).
  * [25] L. M. Alarcón, J. M. Montes de Oca, S. R. Accordino, J. A. Rodriguez Fris and G. A. Appignanesi. Fluid Phase Equilibria 362, 81 (2014).
  * [26] C. A. Menéndez, S. R. Accordino, D. C. Gerbino and G. A. Appignanesi. PLoS One 11 (10), e0165767 (2016).
  * [27] L. M. Alarcón, et al. Eur. Phys. J. E 39, 1-9 (2016).
  * [28] G. Hummer, J. C. Rasaiah and J. P. Noworyta. Nature 414, 188 (2001).
  * [29] N. Giovambattista, P. G. Debenedetti, C. F. Lopez and P. J. Rossky. Proc. Natl. Acad. Sci. U. S. A. 105, 2274 (2008).
  * [30] S. R. Accordino, D. C. Malaspina, J. A. Rodríguez Fris and G. A. Appignanesi. Phys. Rev. Lett. 106, 029801 (2011).
  * [31] J. M. Montes de Oca, F. Sciortino and G. A. Appignanesi. J. Chem. Phys. 152, 244503 (2020).
  * [32] N. A. Loubet, A. R. Verde, J. A. Lockhart, G. A. Appignanesi. J. Chem. Phys. 159, 064512 (2023); doi: 10.1063/5.0159060.
  * [33] J. M. Montes de Oca, S. R. Accordino, A. R. Verde, L. M. Alarcón and G. A. Appignanesi. Phys. Rev. E 99, 062601 (2019).
  * [34] H. J. C. Berendsen, D. van der Spoel and R. van Drunen. Comput. Phys. Commun. 91, 43 (1995).
  * [35] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak. J. Chem. Phys. 81, 3684–3690 (1984).
  * [36] D. A. Case, et al., AMBER9, University of California, SanFrancisco, CA, 2006.