# Characterization of micro-Capsules Deformation in Branching Channels

A. Coclite<EMAIL_ADDRESS>M. D. de Tullio
<EMAIL_ADDRESS>G. Pascazio<EMAIL_ADDRESS>T.
Politi<EMAIL_ADDRESS>Dipartimento di Ingegneria Elettrica e
dell’Informazione, Politecnico di Bari, Via Re David 200 – 70125 Bari, Italy
Dipartimento di Meccanica, Matematica e Management, Politecnico di Bari, Via
Re David 200 – 70125 Bari, Italy

###### Abstract

In this paper, the dynamic of inertial capsules into microfluidic bifurcations
is studied. The fluid evolution is based on the solution of the BGK – lattice
Boltzmann scheme including a forcing term accounting for immersed geometries.
The dynamic-Immersed Boundary forcing strategy is adopted for imposing no-slip
boundary conditions on moving deformable or rigid structures, while, on fixed
immersed geometries the Bouzidi-Firdaouss-Lallemand second-order bounce back
technique is implemented. The proposed computational framework is employed to
detail dynamics and deformation of rigid and deformable capsules traveling
into a branching duct. This journey is characterized in terms of i) the
capsule/bifurcation interaction depending on the sharpness of the branching
channels junction; ii) daughter branches aperture angle; iii) occlusion ratio,
the ratio between capsule size and main channel diameter; iv) flowing capsules
stiffness; v) number of flowing particles.

###### keywords:

dynamic Immersed Boundary , lab on a chip , microfluidics bifurcations , cell
transport , moving least squares

††journal: Applied Mathematics and Computations

## Introduction

As widely known, the micro-circulatory network can be assimilated to a
collection of branching channels in a self-similarity paradigm up to the
micrometric scale [1]. In this context, the transport of a particle dragged by
an external fluid through a branching microchannel is a prototypical problem
that still needs to be thoroughly addressed. The importance of this
fundamental problem resides in the large number of physical phenomena that can
be enlightened depending on the established flow and navigating particles’
characteristics. Such phenomena need to be accurately accounted for describing
the nonlinear behavior of microfluidic devices, such as those used in particle
sorting lab-on-chip[2, 3]. For example, for moderate Reynolds numbers, the
path selected during the journey can be predicted for simple channel
geometries as a function of the particle releasing position and inertia [4, 5,
6]. On the other side, the transport of soft particles in branching ducts is
also considered for the rational design of biomedical devices able to enrich
the concentration of micro-objects with specific characteristics as a result
of the sorting induced by the established flow [7, 8, 9, 10, 11].
Specifically, one can think to focus on the altered mechanical properties and
dynamics due to the altered structure of the cytoskeleton in diseased cells
for the design of such sorting devices [12, 13]. Indeed, such separation
techniques, usually referred as label-free, do not restrict to cells and are
largely used for capsules and particles sorting in microfluidics chip [14, 15,
16, 17]. Recently, the microdevices design community is pushing toward the
idea of mimicking the authentic complexity of nature with in-vitro models for
extravascular mass transport and vascular transport [18, 19, 20]. As a
consequence, the scientific community invested more and more efforts in
reliable computational tools to model the fluid-structure interaction of
micrometric constructs navigating complex geometries mimicking the circulatory
system [21, 22, 23, 24, 25, 26, 27]. All of the considered applications rely
on the specific modeling of the physical phenomena involved and, thus, the
computational tools aiming for the prediction of such phenomena must be
formulated establishing a compromise between accuracy and affordability of
necessary computations.

The present paper aims to discuss the transport and to measure the deformation
of inertial capsules with different stiffness and sizes navigating symmetric
branching vessels as a function of the bifurcation geometry and particles
characteristics. In particular, the authors thoroughly detail the strain
experienced by such capsules underlying the parameters ranges for which large
as well as minimal or negligible deformations are measured. To this scope, a
combined Lattice-Boltzmann / dynamic-Immersed–Boundary (IB) method [28, 29] is
employed due to its ability to predict the transport of massive capsules
within arbitrarily shaped immersed geometries. Specifically, immersed objects
are described by their surface considered as an ensemble of Lagrangian
vertices composing a triangular mesh that moves on the Eulerian grid used to
solve the flow equations. The fluid evolution is obtained through a BGK-
lattice Boltzmann scheme for the Navier-Stokes equation in the incompressible
regime. Dirichlet Boundary Conditions (BC) are enforced on external boundaries
adopting the Zou-He bounce back procedure [30]; while, on immersed objects,
the BCs are imposed in two different fashions depending on whether they move
or not. For fixed obstacles, the Bouzidi-Firdaouss-Lallemand technique is
adopted to bounce back the distribution functions belonging to surrounding
fluid nodes and incoming on solid nodes [31]. On moving geometries, BCs are
enforced by projecting the forcing term computed into the dynamic-IB framework
onto the computational molecula. Structures transport is readily obtained by
solving the equation of motion along with the constitutive model adopted for
particle deformations [29]. The proposed model is implemented to detail the
transport of inertial capsules slightly denser ($\simeq 10\%$) than the
surrounding fluid within a circular cross-section branching channel. This
analysis is carried out systematically by studying the journey of such
capsules in terms of i) capsule/bifurcation interaction and sharpness of the
bifurcation junction and - two different junctions are proposed and analyzed;
ii) symmetric daughter branches aperture angle concerning the main channel
axis $\theta$ (= $\pi/6$, $\pi/3$, and $\pi/2$); iii) occlusion ratio, the
ratio between capsule size $d_{p}$ and main channel diameter $d$ ($d_{p}/d$ =
1/3, 1/2, 2/3, and 5/6); iv) flowing capsules mechanical stiffness, regulated
by the parameter Ca corresponding to the relative effect of hydrodynamics drag
over capsule’s strain resistance - five different values for Ca are considered
(Ca = 0.01, 0.025, 0.05, 0.075, and 0.1); v) number of flowing particles -
single particle as well as multiple particle transport is considered for soft
and rigid capsules.

##  1 Numerical method and computational procedure

In this section a brief description of the dynamic-Immersed–Boundary method
combined with a BGK-Lattice–Boltzmann technique, developed and extensively
validated by the authors is given [28, 29]. The fluid evolution is obtained on
a three-dimensional lattice while the immersed body surface is modeled as a
collection of Lagrangian points subject to a constitutive equation.

###  1.1 Fluid evolution

Figure 1: The D3Q19 lattice molecule. Definition of the D3Q19 computational
molecule with ${\bf e}_{i}$ vectors as defined in Eq.(1) (not in scale).

The incompressible Navier-Stokes equation for the fluid evolution is recovered
with second order accuracy through the BGK-approximation of the Boltzmann
equation on a three-dimensional manifold [32]. Such equation is solved onto
the velocities space on a 3D molecule with 19 reticular velocities namely,
$\\{{\bf e}_{i}({\bf x})\\}_{i=0,\dots,18}$ where ${\bf x}=(x,y,z)\hat{\bf x}$
is the coordinates vector (see Figure.1):

${\bf
e}=\setcounter{MaxMatrixCols}{19}\begin{pmatrix}0&1&-1&0&0&0&0&1&1&1&1&-1&-1&-1&-1&0&0&0&0\\\
0&0&0&1&-1&0&0&1&-1&0&0&1&-1&0&0&1&1&-1&-1\\\
0&0&0&0&0&1&-1&0&0&1&-1&0&0&1&-1&1&-1&1&-1\\\ \end{pmatrix}\,.$ (1)

The solution of the Boltzmann equation is obtained in terms of a vector in the
velocity space composed by 19 discrete functions ${\bf f}=\\{f_{i}({\bf
x})\\}_{i=0,\dots,18}$ evolving on this computational molecula,

${f_{i}({\bf x}+{\bf e}_{i}\Delta t,t+\Delta t)-f_{i}({\bf x},t)=-\frac{\Delta
t}{\tau}[f_{i}({\bf x},t)-f_{i}^{eq}({\bf x},t)]+\Delta t{\cal F}_{i}}\,,$ (2)

where $\tau$ is the relaxation time, $\Delta t$ the time step,
$\\{f^{eq}_{i}({\bf x})\\}_{i=0,\dots,18}$ the components of the equilibrium
discrete functions while ${\cal F}_{i}$ the ones of the discrete forcing term.
The kinematic viscosity $\nu$ is obtained as a function of $\tau$,
$\nu=c_{s}^{2}(\tau-\frac{\Delta t}{2})$ being $c_{s}=1/\sqrt{3}$. Macroscopic
thermodynamics quantities are obtained as statistical moments of ${\bf f}$:

  * 1.

the fluid density $\rho=\sum_{i}f_{i}$;

  * 2.

the momentum $\rho{\bf u}=\sum_{i}f_{i}{\bf e}_{i}+\frac{\Delta t}{2}{\bf
f}_{ib}$;

  * 3.

the thermodynamics pressure $p=c_{s}^{2}\sum_{i}f_{i}=c_{s}^{2}\rho$.

To recover the Navier-Stokes equation, the vector ${\bf f^{eq}}$ is defined as
the expansion of the Maxwell-Boltzmann distribution,

$f_{i}^{eq}({\bf x},t)=\omega_{i}\rho\left[1+\frac{1}{c_{s}^{2}}({\bf
e}_{i}\cdot{\bf u})+\frac{1}{2c_{s}^{4}}({\bf e}_{i}\cdot{\bf
u})^{2}-\frac{1}{2c_{s}^{2}}{\bf u}^{2}\right]\,,$ (3)

with $\omega_{0}=1/3$, $\omega_{1-6}=1/18$, and $\omega_{7-18}=1/36$. In the
same fashion, the expansion of ${\cal F}_{i}$ reads [33, 34]:

${{\cal F}_{i}=\Bigl{(}1-\frac{1}{2\,\tau}\Bigr{)}\omega_{i}\Bigl{[}\frac{{\bf
e_{i}}-{\bf u}}{c_{s}^{2}}+\frac{{\bf e_{i}}\cdot{\bf u}}{c_{s}^{4}}{\bf
e_{i}}\Bigr{]}\cdot{\bf f}_{ib}}\,$ (4)

where ${\bf f}_{ib}({\bf x},t)$ is the force term due to the presence of
immersed structures [35, 28].

###  1.2 Dirichlet conditions for external boundaries

Dirichlet Boundary Conditions (DBC) are enforced by adopting the bounce back
procedure by Zou and He [30] to enforce velocity- or pressure-based boundary
conditions. As widely known Eq.(2) is solved numerically in two steps, the
collision, and the streaming phases. However, for boundaries, the streaming
phase cannot be completed since it is not possible to stream in and out the
distribution functions for the boundary nodes. Let us consider, for example,
the north wall ($z=z_{max}$). All of the distribution functions entering the
wall are known since they derive from the post-collision streaming phase of
the nodes belonging from the plane $z=z_{maz}-\Delta x$, namely
$f_{5},f_{9},\,f_{11},\,f_{16}$ and $f_{18}$ (see Figure.1). On the contrary,
all of the distribution functions that should be streamed back into the wall
nodes from the plane $z=z_{max}+\Delta x$, namely
$f_{6},f_{10},\,f_{12},\,f_{15}$ and $f_{17}$ (see Figure.1) are unknown and
need to be somehow imposed or computed. Zou and He proposed to enforce
conservative boundary conditions through a sophisticated bounce-back procedure
in which the unknown distribution functions are computed also by imposing the
value of three macroscopic quantities, namely, the three components of the
velocity field for velocity-based DBCs or the pressure and the two non-
orthogonal components (with respect to the boundary normal) of the velocity
for pressure-based DBCs. Therefore, for the north wall the velocity-based DBCs
read:

${\begin{cases}\rho(x_{w},y_{w},z_{w})=(u_{z,w}+1)^{-1}\Bigl{[}f_{0}+f_{1}+f_{2}+f_{3}+f_{4}+f_{7}+f_{8}+f_{13}+f_{14}+2\Bigl{(}f_{5}+f_{9}+f_{11}+f_{16}+f_{18}\Bigr{)}\Bigr{]}\\\
f_{6}=f_{5}-\frac{1}{3}\rho(x_{w},y_{w},z_{w})u_{z,w}\\\
f_{15}=f_{16}+\frac{1}{6}(-u_{z,w}+u_{x,w})-N_{z,x}\\\
f_{10}=f_{9}+\frac{1}{6}(-u_{z,w}-u_{x,w})+N_{z,x}\\\
f_{17}=f_{18}+\frac{1}{6}(-u_{z,w}+u_{y,w})-N_{z,y}\\\
f_{12}=f_{11}+\frac{1}{6}(-u_{z,w}-u_{y,w})+N_{z,y}\\\ \end{cases}}\,,$ (5)

where $u_{x,w}$, $u_{y,w}$, and $u_{z,w}$ are the values of the velocity
components imposed at the wall nodes and

$\displaystyle N_{z,x}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigl{(}f_{1}+f_{7}+f_{13}-(f_{2}+f_{14}+f_{8})\Bigr{)}-\frac{1}{3}\rho(x_{w},y_{w},z_{w})u_{x,w}\,,$
(6) $\displaystyle N_{z,y}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\Bigl{(}f_{3}+f_{7}+f_{14}-(f_{4}+f_{13}+f_{8})\Bigr{)}-\frac{1}{3}\rho(x_{w},y_{w},z_{w})u_{y,w}\,.$
(7)

###  1.3 No-slip conditions for moving immersed boundaries

In the Immersed-Boundary technique [36], the surface of an arbitrarily shaped
body immersed in the flow is modeled as an ensemble of Lagrangian vertices
moving on the withstanding Eulerian grid on which the fluid solution is
computed (Figure.2a). Here, such surface is discretized utilizing a triangular
mesh and the centroids of each element are called “Lagrangian markers”. A
convenient transfer function is built on a moving least squares reconstruction
for exchanging needed information between the two meshes [28, 29, 37].

Figure 2: Immersed-boundary technique schematic. a. Spatial approximation of
an immersed body obtained with a triangular mesh. Each element centroid
corresponds to a Lagrangian marker. b. Eulerian points (per plane) on which
the dynamic-IB forcing procedure is obtained through the support (blue square)
centered on the selected Lagrangian marker (green circle). c. Eulerian points
(per plane) used for hydrodynamics stresses computations taken at 1.2 $\Delta
x$ along the inward and outward normal directions of each element.

Specifically, for each triangular element centroid, 9 Eulerian points per
plane are identified as all of the Eulerian points falling into the support
domain by side equal to $2r_{w}=2.6\,\Delta x$ and centered on the selected
centroid (indexed here with $l$); $\Delta x$ being the Eulerian grid spatial
resolution (see Figure.2b). Given the fluid solution at time step $n$, the
velocity the marker $l$ is approximated by,

${U_{i}({\bf x})=\sum_{k=1}^{27}\,\phi^{l}_{k}({\bf x})u^{k}_{i}}\,,$ (8)

being $u^{k}_{i}$ the component $i$ of the fluid velocity field ${\bf u}$
taken at the Eulerian point $k$ associated to the marker $l$. $\phi$ is the
operator given by the minimization with respect to ${\bf a(x)}$ of:

${J=\sum_{k=1}^{27}W({\bf x}-{\bf x}^{k}})[{\bf p}^{T}({\bf x}^{k}){\bf
a(x)}-u^{k}_{i}]^{2}\,.$ (9)

In the latter, ${\bf x}$ is the coordinate vector of the marker, ${\bf p}^{T}$
is a linear basis vector, ${\bf a(x)}$ a vector of coefficients such that
$\sum_{k=1}^{27}\,\phi^{l}_{k}({\bf x})u^{k}_{i}\,={\bf p}^{T}({\bf
x}^{k}){\bf a(x)}$. $W({\bf x}-{\bf x}^{k})$ is the weight function:

$W({\bf x}-{\bf x}^{k})=\begin{cases}e^{-(r_{k}/\alpha)^{2}}\,&\mbox{if
}r_{k}\leq 1\\\ 0\,&\mbox{if }r_{k}>1\\\ \end{cases}\,,$ (10)

where, $\alpha=0.3$ and $r_{k}$ is the distance between the marker and point
$k$ normalized with respect to the support domain size length. For each
Lagrangian marker, the volume force ${\bf F}_{l}$ required to enforce no-slip
BCs is computed, and then transferred to the Eulerian points falling into the
conveniently defined support domain. Once the velocity at the $l$-th
Lagrangian marker is known from Eq. 8, the desired velocity at that point,
${\bf U}_{b}^{l}({\bf x})$, is enforced by the volume force:

${{\bf F}_{l}{\bf(x)}=\frac{{\bf U}_{b}^{l}{\bf(x)}-{\bf U(x)}}{\Delta t}}\,,$
(11)

that is then used to compute the force for the point $k$,

${{\bf f}_{ib}}^{k}=\sum_{l}c_{l}\phi^{l}_{k}{\bf F}_{l}\,.$ (12)

Note that, in the above equation the sum runs over all markers having the
point $k$ point in their support domain. $c_{l}$ is a coefficient obtained by
imposing the conservation of the total force exerted by the fluid [38] between
the two grids.

###  1.4 Immersed-Structure Transport

Immersed bodies motion due to the interaction established with the fluid is
modeled by adopting the dynamics-IB technique described in [28, 29]. For the
case of soft objects, the equation of motion of each Lagrangian point is
solved considering both constitutive model response and hydrodynamics forces.
The case of rigid bodies is treated separately by integrating Newtown’s second
law and the moment equation for each body.

Deforming bodies constitutive model. Deforming bodies internal forces are
computed as elastic strain response, bending resistance, and enclosed volume
conservation [29]. The stretching elastic potential acting on the two vertices
sharing the $l$-th edge is given by

${V_{l}^{s}=\frac{1}{2}k_{s,l}(l_{l}-l_{l,0})^{2}}\,,$ (13)

being $k_{s,l}$ [N/m] the elastic modulus for the edge $l$, $l_{l}$ its length
in the deformed configuration and $l_{l,0}$ its length at the resting. The
elastic forces for vertices 1 and 2 connected by the edge l are:

${\begin{cases}{\bf F}_{1}^{int,s}=-k_{s,l}(l-l_{0})\frac{{\bf
r}_{1,2}}{l}\,\\\ {\bf F}_{2}^{int,s}=-k_{s,l}(l-l_{0})\frac{{\bf
r}_{2,1}}{l}\,\\\ \end{cases}}\,,$ (14)

where ${\bf r}_{i,j}={\bf r}_{i}-{\bf r}_{j}$ and ${\bf r}_{i}$ is the ray
vector connecting i to j. The elastic modulus of the immersed surface is
obtained trough the nondimensional group measuring the ratio between viscous
drag forces versus surface tension forces, the capillary number
Ca$=\frac{\rho\nu u_{ref}}{k_{s,0}}$. Then, it is spread to all edges
considering the fraction of the total surface assigned to each edge, namely
the sum of the area of the two elements sharing the selected edge $A^{l}_{i}$
: $k_{s,l}=\frac{k_{s,0}T\sum_{i}A^{l}_{i}}{l^{2}}$ [39]. Indeed, $\rho$ and
$\nu$ are the fluid density and kinematic viscosity; $u_{ref}$ the reference
velocity, $k_{s,0}$ the surface elastic modulus, and $T$ the surface
thickness.

Off-plane deformations are constrained by a bending potential depending on the
difference between the local curvature in the resting and the deformed
configurations, $c_{v,0}$ and $c_{v}$, respectively,

${V_{v}^{b}=\frac{1}{2}k_{b}(c_{v}-c_{v,0})^{2}}\,,$ (15)

being $k_{b}$ the bending modulus. $c_{v}$ is computed through the angle
defined by the unit normal vectors of two adjacent elements, ${\bf n}_{1}$ and
${\bf n}_{2}$, ($\theta-\theta_{0}$), with $\theta_{0}$ the angle in the
resting configuration. $k_{b}$ is obtained as a function of $k_{s,l}$ through
the nondimensional group $Eb=\frac{k_{b}}{k_{s,l}r^{2}}$. [40] The resulting
force on the involved four nodes reads:

${\begin{cases}{\bf F}^{int,b}_{1}=\beta_{b}\bigl{[}b_{11}({\bf
n}_{1}\times{\bf r}_{32})+b_{12}({\bf n}_{2}\times{\bf r}_{32})\bigr{]}\,\\\
{\bf F}^{int,b}_{2}=\beta_{b}\bigl{[}b_{11}({\bf n}_{1}\times{\bf
r}_{13})+b_{12}({\bf n}_{1}\times{\bf r}_{34}+{\bf n}_{2}\times{\bf
r}_{13})+b_{22}({\bf n}_{2}\times{\bf r}_{34})\bigr{]}\,\\\ {\bf
F}^{int,b}_{3}=\beta_{b}\bigl{[}b_{11}({\bf n}_{1}\times{\bf
r}_{21})+b_{12}({\bf n}_{1}\times{\bf r}_{42}+{\bf n}_{2}\times{\bf
r}_{21})+b_{22}({\bf n}_{2}\times{\bf r}_{42})\bigr{]}\,\\\ {\bf
F}^{int,b}_{4}=\beta_{b}\bigl{[}b_{11}({\bf n}_{1}\times{\bf
r}_{23})+b_{22}({\bf n}_{2}\times{\bf r}_{23})\bigr{]}\,\\\ \end{cases}}\,,$
(16)

with $b_{ij}=\frac{{\bf n}_{i}\cdot{\bf n}_{j}}{|{\bf n}_{i}||{\bf n}_{j}|}$
and
$\beta_{b}=k_{b}\frac{sin(\theta)cos(\theta_{0})-cos(\theta)sin(\theta_{0})}{\sqrt{1-cos^{2}(\theta)}}$.

The incompressibility (conservation of the total volume) of the membrane is
enforced through a penalty model:

${{\bf F}_{i}^{int,v}=-k_{v}\left(1-\frac{V}{V_{0}}\right){\bf
n}_{i}A_{i}}\,,$ (17)

$k_{v}$ an incompressibility constant; $V$ and $V_{0}$ the current and initial
volume enclosed by the structure; $A_{i}$ the area associated with the element
$i$.

Hydrodynamics stresses. The hydrodynamic force acting on an element is
composed by two terms, ${\bf F}_{i}^{ext,p}(t)$ and ${\bf
F}_{i}^{ext,\tau}(t)$:

${{\bf F}_{i}^{ext,p}(t)=-(p^{+}_{i}-p^{-}_{i})\,{\bf n}^{+}_{i}\,A_{i}}\,,\\\
$ (18) ${{\bf
F}_{i}^{ext,\tau}(t)=(\bar{\bf\tau}^{+}_{i}-\bar{\bf\tau}^{-}_{i})\cdot{\bf
n}^{+}_{i}\,A_{i}}\,.$ (19)

$p^{+}_{i}$, $p^{-}_{i}$ and $\bar{\bf\tau}^{+}_{i}$, $\bar{\bf\tau}^{-}_{i}$
are the pressure and the viscous stress tensor taken on the marker associated
to element $i$. Such, quantities accounts for both external (+) and internal
(-) flow fields. ${\bf n^{+}}_{i}$ is the normal outward unit vector. Pressure
and velocity derivatives in Eq.s (18) and .(19) are computed at a distance 1.2
$\Delta x$ along the outward (+) and inward (-) normal directions (see
Figure.2c) [37].

Bodies interaction and collision treatment. The interaction between two or
more immersed bodies is treated by imposing a purely repulsive potential
centered in each Lagrangian marker of each immersed object. Such repulsive
force is such that the minimum distance allowed between two Lagrangian markers
belonging to two different immersed structures is $\Delta x$. Specifically,
the force exerted by the $i$-th Lagrangian marker of a particle at a distance
$d_{i,j}$ from the $j$-th marker of a second particle is chosen to be
proportional to the module of the characteristic velocity of the flow
$u_{max}$:

${{\bf F}_{i}^{ext,c}=\frac{u_{max}\Delta t}{8\sqrt{2}}\sqrt{\frac{\Delta
x}{d_{i,j}^{5}}}{\bf n}^{-}_{i}}\,.$ (20)

Soft structure motion. The position of each element centroid is updated in
time by integrating the equation of motion:

${m_{v}\dot{{\bf u}}_{v}={\bf F}_{v}^{int,s}(t)+{\bf F}_{v}^{int,b}(t)+{\bf
F}_{v}^{int,v}(t)+{\bf F}_{v}^{ext,c}(t)+{\bf F}_{v}^{ext,p}(t)+{\bf
F}_{v}^{ext,\tau}(t)}\,,$ (21)

where $m_{v}$ is the mass associated to the point $v$ with the assumption that
the structure mass is uniformly distributed over the vertices; $\dot{{\bf
u}}_{v}$ corresponds to the Lagrangian point acceleration. Such equation is
integrated in time adopting the Verlet algorithm with the zero-step tentative
velocity ${\bf u}_{v,0}(t+\Delta t)$ obtained by interpolating the fluid
velocity from the Eulerian points falling in the support domain centered in
$v$,

${{\bf x}_{v}(t+\Delta t)={\bf x}_{v}(t)+{\bf u}_{v,0}(t+\Delta t)\Delta
t+\frac{1}{2}\frac{{\bf F}_{v}^{tot}(t)}{m_{v}}\Delta t^{2}+O(\Delta
t^{3})}\,;$ (22)

then the vertex velocity is computed as:

${{\bf u}_{v}(t+\Delta t)=\frac{\frac{3}{2}{\bf x}_{v}(t+\Delta t)-2{\bf
x}_{v}(t)+\frac{1}{2}{\bf x}_{v}(t-\Delta t)}{\Delta t}+O(\Delta t^{2})}\,,$
(23)

and ${\bf x}_{v}(t+\Delta t)$ is finally corrected:.

${{\bf x}_{v}(t+\Delta t)={\bf x}_{v}(t)+{\bf u}_{v}(t+\Delta t)\Delta
t+\frac{1}{2}\frac{{\bf F}_{v}^{tot}(t)}{m_{v}}\Delta t^{2}+O(\Delta
t^{3})}\,.$ (24)

Rigid structures transport. The transport of rigid structures is obtained by
summing all force contributions over the structure surface (Eq.s.(18), .(19)
and .(20)) and then computing the linear and angular velocities and positions.
The Newton-Euler system reads,

$\begin{cases}{\bf F}^{tot}(t)=m\,\dot{{\bf u}}(t)\\\ {\bf M}^{tot}(t)={\bf
I}\,\dot{{\bf\omega}}(t)\end{cases}\,.$ (25)

${\bf F}^{tot}(t)$ and ${\bf M}^{tot}(t)$ are the total force and moment; $m$
is the total structure mass; ${\bf I}$ the inertia tensor. ${\bf u}(t)$ and
${\bf\omega}(t)$ are derived as:

$\begin{cases}{\bf u}(t)=\frac{2}{3}(2{\bf u}(t-\Delta t)-\frac{1}{2}{\bf
u}(t-2\Delta t)+\dot{{\bf u}}(t)\Delta t)+O(\Delta t^{2})\\\
{\bf\omega}(t)=\frac{2}{3}(2{\bf\omega}(t-\Delta
t)-\frac{1}{2}{\bf\omega}(t-2\Delta t)+\dot{{\bf\omega}}(t)\Delta t)+O(\Delta
t^{2})\end{cases}\,.$ (26)

###  1.5 No-Slip boundary conditions for stationary immersed boundaries

Figure 3: Bouzidi-Firdaous-Lallemand interpolation technique. Definition of
solid and boundary nodes involved into the enforcement of the boundary
conditions along with the distance q used for the interpolation of the
bounced-back distribution functions.

To limit the computational burden and save computational time, boundary
conditions for stationary bodies and moving geometries are treated with two
different techniques. Specifically, the second-order Bouzidi-Firdaouss-
Lallemand technique is adopted to bounce-back distribution functions belonging
from surrounding computational nodes incoming against the computational nodes
occupied by the solid obstacle [31]. In the same fashion used for moving
geometries, in this context, an immersed body corresponds to a collection of
triangular elements. All of the computational nodes are tagged as “fluid”,
“solid” or “boundary”. Eulerian points corresponding to computational cells in
which a Lagrangian marker falls are considered “solid nodes”. Then, a node
into the computational domain not tagged as solid that is at a distance ${\bf
e}_{i}\Delta t,\,i=0,...,18$ from a solid node is tagged as “boundary” node,
being ${\bf e}_{i}$ the 19 reticular directions, all other nodes are
considered as “fluid nodes” (see Figure.3). Note that, within this
parametrization, boundary nodes play a major role being the nodes on which the
no-slip boundary conditions are enforced. Once the computational domain is
completely tagged, for all of the boundary nodes are identified the reticular
directions along which at least a solid point is found and unity value is
assigned to the scalar field $a_{i}({\bf x})$ otherwise null (see Figure.3).
Moreover, for such directions, the normalized distance $q_{i}({\bf x})$ is
computed as the ratio between the distance from the corresponding Lagrangian
marker (with coordinates $({\bf x}_{L})$) and the boundary node (with
coordinates $({\bf x}_{b})$) over the distance between the boundary node and
the solid node encountered along the $i$-th reticular direction (with
coordinates $({\bf x}_{s})$):

${q_{i}({\bf x})=\frac{\sqrt{({\bf x}_{b}-{\bf x}_{L})^{2}}}{\sqrt{({\bf
x}_{b}-{\bf x}_{s})^{2}}}}\,.$ (27)

Finally, for all of the boundary nodes the incoming post-collision
distribution functions are bounced-back following the rule:

${\begin{cases}f_{\overline{i}}({\bf x})=q_{i}({\bf x})(2q_{i}({\bf
x})+1)f_{i}({\bf x})+(1-4q_{i}^{2}({\bf x}))f_{i}({\bf x}-{\bf e}(i)\Delta
t)-q_{i}({\bf x})(1-2q_{i}({\bf x}))f_{i}({\bf x}-2{\bf e}(i)\Delta t)\\\
$\hfill if\quad$q_{i}({\bf x})<0.5\\\ f_{\overline{i}}({\bf
x})=\frac{1}{q_{i}({\bf x})(2q_{i}({\bf x})+1)}f_{i}({\bf
x})+\frac{2q_{i}({\bf x})-1}{q_{i}({\bf x})}f_{i}({\bf x})+\frac{1-2q_{i}({\bf
x})}{1+2q_{i}({\bf x})}f_{i}({\bf x}-{e}_{i}\Delta t)$\hfill
if\quad$q_{i}({\bf x})\geq 0.5\par\end{cases}}\,,$ (28)

where $\overline{i}$ indicates the opposite reticular direction with respect
to $i$.

##  2 Results and discussion

###  2.1 Computational set-up and boundary conditions

The deformation of a capsule navigating a microfluidic bifurcation is studied
here in terms of the capsule mechanical stiffness Ca, the occlusion ratio
$d_{p}/d$ (ratio between the particle $d_{p}$ and the main channel $d$
diameters), the bifurcation aperture angle $\theta$ as well as its junction
geometry, and the number of flowing capsules. In particular, this systematical
analysis is composed of two numerical experiments campaigns: i) the transport
of a single capsule and ii) the transport of a train of three capsules into a
microfluidic bifurcation. Firstly, the authors assess the transport of an
isolated capsule by measuring the relative variation of the capsule surface
area $\delta S$ while navigating the microfluidic chip for all of the possible
combinations of the geometrical and mechanical parameters involved. The
ability of such isolated capsule of selecting a path toward one of the two
branches is not considered in this work. As a matter of fact, due to the
length of the main channel, it is reasonable to suppose that such capsule
would reach the channel midway before entering the junction region. With this
rationale, the releasing lateral position is kept constant as 0.5 d. Note
that, for small Reynolds number, hydrodynamics do not prescribe any reason to
choose one of the two daughter branches. The role of the capsule’s inertia is
non-negligible only during the interaction with the vasculature. Due to these
rationales, the authors do not discuss the path-selection problem since it is
random in principle. Indeed, the randomness of nature losts in computations
and the selection of one of the two branches may be due to several operating
factors like the indexing of Lagrangian mesh triangles in exerted stresses
evaluations or Eulerian grid resolution. Then, for a specific occlusion ratio,
$d_{p}/d=2/3$, the journey of a train of three capsules as a function of Ca
and $\theta$ is investigated and characterized in terms of the path traveled
and relative stretching. Lastly, this journey is assessed when transporting
smaller rigid particles within the train of capsules and analyzing their
mutual interaction.

A circular cross-section channel with diameter d = 30 $\Delta x$ branches in
two daughter ducts with diameter 0.7 d so as to maintain constant the total
cross-section (see Figure.4a). The main channel length is 6.7 d while the two
branches are 5 d long each. The semi-aperture angle $\theta$ is computed as
the angle between the main channel axis and the axis of one of the two
branches and corresponds to $\theta=$ $\pi/6$, $\pi/3$, and $\pi/2$. Moreover,
two different geometries for the junction of the three channels are
considered, a sharp junction and a smoother joint with a trapezoidal vertical
section, reported in Figure.4. An initially spherical capsule is placed at
1.16 d from the inlet section (see Figure.4a); 4 different values of the
particle diameter are considered, namely, $d_{p}=\frac{1}{3}d$,
$\frac{1}{2}d$, $\frac{2}{3}d$, and $\frac{5}{6}d$. The number of triangular
elements $n_{L}$ discretizing the capsule is 325, 485, 645 and 785 for
$d_{p}/d$ = 1/3, 1/2, 2/3, and 5/6, respectively. As detailed in the
computational method section, the mechanical stiffness of such object is
regulated by two parameters, the capillary number Ca and the bending
resistance Eb. In this work Ca = 0.01, 0.025, 0.05, 0.075 and 0.1 while the
bending modulus equals Eb = 0.001. The capsules mass is computed through the
ratio between solid and fluid densities, fixed here as
$\rho_{s}/\rho_{f}=1.1$.

Figure 4: Schematic of the microfluidic bifurcation. a. Sketch of the
microfluidic chip along with characteristic dimensions and lengths when
considering the sharp junction between the main and the two daughter-branches.
(b, c, d) Sketch of the microfluidic bifurcations with the smooth junction for
$\theta=\pi/6$ (b), $\theta=\pi/3$ (c), and $\theta=\pi/2$ (d).

At the inlet section, the velocity-based Dirichlet conditions are imposed and
the Poiseuille parabolic velocity profile with peak velocity $u_{max}$ is
adopted, thus, mimicking the velocity profile of an already established flow.
At the two outlet sections the Dirichlet pressure-based outlet conditions with
$p_{out}=p_{ref}$ are imposed. $u_{max}$ is regulated by the Reynolds number
of the flow defined as Re = $\frac{u_{max}d}{\nu}$ = 0.01. The kinematic
viscosity $\nu$ depends on the relaxation time as
$\nu=c_{s}^{2}(\tau-\frac{\Delta t}{2})$ being $\tau=0.75$ and $\Delta t=1$
[29]. Notably, Re = 0.01 correspond to a peak velocity of $u_{max}$ = 0.001
m/s for a 10 $\mu$m particle transported in plain water. The flow field is
initialized with the solution of the stationary flow obtained when no
particles or capsule is transported into the bifurcation (see Supp.Figure.1
and .2).

###  2.2 Transport of an elastic capsule in a microfluidic bifurcation

The role of the junction geometry. The transport of the capsule with $d_{p}/d$
= 2/3 and the analysis of its deformation when crossing the microfluidic chip
is firstly considered to compare the abilities of the two junctions proposed
in Figure.4a and .4b. Specifically, the ability of the two junctions in
facilitating the journey of the microcapsule is assessed considering the
probability of large deformations induced by the specific geometry. In this
work, a link connecting two vertices of the Lagrangian mesh composing the
capsule is considered damaged or broken when the actual length doubles its
value in the stress-free configuration (at this occurrence the computation is
stopped). The deformation of a single capsule transported within the
bifurcation as a function of the mechanical stiffness for $\theta=\pi/6$ is
analyzed. The configurations acquired during the transport of the capsule for
Ca = 0.01, 0.05 and 0.1 are documented for the sharp and the smooth junction
in Figure.5a and in 5b, respectively.

Figure 5: Single capsule crossing the bifurcation for $\theta=\pi/6$ as a
function of Ca. Configurations of a single capsule navigating a bifurcation
with the sharp (a.) and the smooth junction (b.) for Ca = 0.01 (left column),
Ca = 0.05 (central column), and Ca = 0.10 (right column). The configurations
correspond to $tu_{max}/d$ from 0 to 40 with a step of 3.5. Close-ups refer to
impinging instant depending, indeed, by the capsule stiffness.

The capsule is released for $tu_{max}/d=0$ and navigates the main channel
while deforming under the established flow. Then, entering the junction region
the capsule impinges on the upper wall and, depending on the geometry, it may
be stretched and damaged. When capsules hit the cusp in the sharp junction
(see close-up in Figure.4a) they squeeze while adhering to the cusp. At this
stage, the capsule peripheries are equally distributed in both daughter
branches (as documented in the close-ups of Figure.5a). Immediately after this
phase, the capsule randomly chooses a path to follow selecting one of the two
branches and crawling the cusp edge eventually navigating the selected duct.
However, for Ca = 0.1 the capsule results damaged during the impingement with
the cusp edge; in fact, it deforms more and more responding to the symmetry of
the flow, which never pushes the particle in one of the two branches.

Figure 6: Trajectory and relative stretching of the single capsule crossing
the bifurcation with $\theta=\pi/6$ as a function of the mechanical stiffness.
(a, b) Path traveled by a capsule into the microfluidic bifurcation for
$\theta=\pi/6$ and Ca = 0.01, 0.025, 0.05, 0.075, and 0.1 with the sharp (a)
and the new (b) junction. Distribution of $\delta S$ over time obtained with
the sharp (c) and the new (d) junction. Distribution of $\delta S$ over the
vertical coordinate $z/d$ for the sharp (e) and the new (f) junction. The
dashed lines correspond to the vertical coordinate of the junction edge.

This picture results quite different when considering the smoother joint
presented in Figure.4b. Indeed, by adopting a junction with a trapezoidal
section, the capsule will impinge and adhere to a flat wall. Then, again, it
deforms responding to the symmetry of the flow, and crawls over this flat wall
until it is dragged into one of the two smaller vessels (see Figure.5b).
Interestingly, with this smooth junction, the capsule can overcome the
bifurcation and travel one of the two branches for all of the investigated
values of Ca. However, at the impinging instant, all of the capsules
transported within the smooth junction present smaller deformation, as
reported in Figure.5b. The trajectory of the transported capsules as a
function of the stiffness is plotted Figure.6 along with the distribution of
the relative stretching $\delta S$. The relative stretching is defined as the
ratio $\frac{S(t)-S_{0}}{S_{0}}$ where $S(t)$ and $S_{0}$ are the current and
the initial capsule surface area. As discussed, the capsules stiffness does
not seem to play any role in the choice of which branch will be traveled after
overcoming the junction. Moreover, the distribution in time of the relative
stretching demonstrates how for a fixed capillary number the geometry of the
junction may induce larger or smaller deformation along with the occurrence of
the capsule membrane degradation. For Ca = 0.05, for example, the peak value
of $\delta S$ is 0.161 (for $tu_{max/d}$ = 14.28) for the sharp junction and
0.058 (for $tu_{max/d}$ = 15.70) for the smooth junction; for Ca = 0.075, the
capsule results damaged when using the sharp junction while with the smooth
junction the peak value of $\delta S$ registered corresponds to 0.22 at
$tu_{max}/d$ = 19.04 (see Figure.6c and .6d). The time instant in which the
peak value of $\delta S$ is registered depends, indeed, on the capsule’s
mechanical stiffness. In fact, on one side, softer capsules spend a part of
the momentum transferred by the fluid for deforming their surface; on the
other side, stiffer capsules are gradually deformed by the flow and spend more
time acquiring the needed deformation to navigate a daughter vessel. On the
contrary, the position along the vertical coordinate in which the most of the
deformation is performed, so that, in which the peak of $\delta S$ is
registered, depends on the peculiar flow patterns obtained and not only on Ca.
However, when considering the sharp junction, the position in which those
peaks are registered varies with Ca, as documented in Figure.6e. While, this
phenomenon becomes negligible whit the smooth junction (see Figure.6f). The
peak values registered when using the sharp junction correspond to 6.22, 6.30,
and 6.39 $z/d$ for Ca = 0.01, 0.025, and 0.05, respectively. While such peaks
are measured at 6.70, 6.72, 6.72, 6.72, and 6.74 $z/d$ for Ca = 0.01, 0.025,
0.05, 0.075, and 0.1 whit the smooth junction.

Figure 7: Flow patterns induced by the presence of the transported capsule
during the impinging with the junction edge for $\theta=\pi/6$ and Ca = 0.01,
0.05 obtained with both junctions.

Such difference is ascribable to the peculiar flow patterns obtained with the
two. The smooth junction presents a flat edge close to which the fluid is
resting like in a sort of stagnation region (see no-cell contour plots in
Figure.7). This stagnation region slows down the capsules letting them to
acquire the needed deformation to navigate the secondary branches. Indeed,
also the fluid flow is perturbed by the capsule presence in two different
fashions when considering the two junctions. Such perturbations are measured
here by seeing at the vertical component of the velocity field $u_{z}/u_{max}$
as well as the relative pressure $(p-p_{0})/p_{0}$, being $p_{0}$ the outflow
pressure (see Figure.7). When the capsule impinges on the cusp edge of the
sharp junction it squeezes adhering to it and implying a partial occlusion of
both branches whose extension depends on Ca. Softer capsules, indeed, occupy
larger regions of both branches during this phase (see the plots for $z/d=6.5$
obtained with the sharp junction in Figure.7). Interestingly, such partial
occlusion is negligible when using the smooth junction and the capsules are
able to navigate the vasculature without spreading into both branches and
perturbing the flow in both secondary channels (see the lower panel in
Figure.7).

Figure 8: Stretching of a single capsule transported in the bifurcation for
$\theta=\pi/6$ as a function of $d_{p}/d$ and Ca for both junctions.

The ability of the two junctions in facilitating the journey of a microcapsule
is further analyzed in terms of the occlusion ratio $d_{p}/d$. Specifically,
the particle diameter is varied from 1/3 d to 5/6 d with a 1/6 d step. The
distributions of $\delta S$ over time are documented in Figure.8 for $d_{p}/d$
= 1/3, 1/2, and 5/6. As expected, for $d_{p}/d$ = 1/3 the capsule experience
almost negligible deformations (smaller than 0.05) for all of the investigated
capillary numbers. However, the sharp junction still implies much stronger
deformations with respect to the new one. For $d_{p}/d$ = 1/2 deformations are
valuable and the peak values registered are almost doubled when considering
the sharp over the smooth junction. Specifically, the peak value of $\delta S$
corresponds to 0.106, 0.154, and 0.235 (taken at $tu_{max/d}$ = 13.33, 13.57,
and 14.28) for Ca = 0.05, 0.075 and 0.1, respectively, when using the sharp
joint; while, it equals 0.025, 0.054, and 0.086 (taken at $tu_{max/d}$ =
16.03, 16.67, and 17.14) for Ca = 0.05, 0.075 and 0.1, respectively, when
using the smooth junction (see plots in the middle row of Figure.8). Notably,
for $d_{p}/d$ = 1/2 no capsule results damaged during the journey. On the
contrary, for the largest occlusion ratio the capsule results damaged for Ca =
0.05, 0.075 and 0.1 for the sharp junction and for Ca = 0.075 and 0.1 for the
new joint. Indeed, the measured peak values of $\delta S$ increase with
$d_{p}/d$ due to needed squeezing for passing the bifurcation. From this
comparison, we conclude that smoother geometries for the junctions should be
preferred in the design phase of a microfluidic chip against a net joining of
the two daughter branches onto the main channel and from now on we will
further analyze the deformation of capsule only considering the smooth
junction.

The role of the bifurcation aperture angle. The deformation of a capsule
navigating a microfluidic bifurcation is further assessed when varying the
semi-aperture angle of the two branches with respect to the main channel (see
Figure.4c and .4d). For $\theta=\pi/3$ and $\pi/2$ the deformation experienced
by the capsule at the impinging instant is smaller than that corresponding to
$\theta=\pi/6$ (see Figure.9a and .9b).

Figure 9: Single capsule transported in the bifurcation for $\theta=\pi/3$ and
$\pi/2$ (a, b) Configuration of the capsule at the impinging instant for Ca =
0.01 (left column), 0.05 (central column), and 0.1 (right column) for
$\theta=\pi/3$ (a) and $\pi/2$ (b). (c, d) Trajectory and relative stretching
of the capsule navigating the bifurcation for $\theta=\pi/3$ (c) and $\pi/2$
(d) as a function of Ca.

Indeed, for $\theta=\pi/6$ the flow velocity has a large component along z in
the junction region dramatically pushing the capsule against the upper wall of
the junction (see Supp.Figure.1); however, for larger $\theta$ the geometry of
the bifurcation implies an abrupt change in the flow direction that results in
a small ($\theta=\pi/3$) or null ($\theta=\pi/2$) component along z into the
junction region and the capsule is essentially dragged by inertia against the
junction boundary (see Supp.Figure.1). The deformed configuration of the
capsule as a function of Ca at the impinging instant is documented in
Figure.9a and 9b for $\theta=\pi/3$ and $\pi/2$, respectively. Note that, for
Ca = 0.01, the capsule does not touch the junction wall for both values of
$\theta$. Moreover, also in these two cases the path of such capsules do not
seem to be influenced by Ca, as demonstrated in Figure.9c and .9d. The
distribution of the relative surface stretching as a function of the stiffness
is reported in Figure.9. Quantitatively, the peak values of $\delta S$
registered for Ca = 0.075 are 0.105 (at $tu_{max}/d$ = 14.3) and 0.101 (at
$tu_{max}/d$ = 16.2) for $\theta=\pi/3$ and $\pi/2$, respectively,
corresponding in both cases at about an half of the peak values obtained for
$\theta=\pi/6$.

|  |  $\theta=\pi/6$ Sharp Jun. |  $\theta=\pi/6$ New Jun. | $\theta=\pi/3$ | $\theta=\pi/2$
---|---|---|---|---|---
| Ca = 0.01 | $\leq 0.002$ | $\leq 0.002$ | $\leq 0.002$ | $\leq 0.002$
| Ca = 0.025 | $\leq 0.002$ | $\leq 0.002$ | $\leq 0.002$ | $\leq 0.002$
$d_{p}/d$ = 1/3 | Ca = 0.05 | 0.016 | $\leq 0.002$ | $\leq 0.002$ | $\leq 0.002$
| Ca = 0.075 | 0.035 | $\leq 0.002$ | $\leq 0.002$ | $\leq 0.002$
| Ca = 0.1 | 0.047 | $\leq 0.002$ | $\leq 0.002$ | $\leq 0.002$
| Ca = 0.01 | 0.017 | 0.006 | $\leq 0.002$ | $\leq 0.002$
| Ca = 0.025 | 0.042 | 0.009 | 0.007 | 0.005
$d_{p}/d$ = 1/2 | Ca = 0.05 | 0.106 | 0.025 | 0.021 | 0.011
| Ca = 0.075 | 0.154 | 0.054 | 0.023 | 0.017
| Ca = 0.1 | 0.235 | 0.086 | 0.028 | 0.024
| Ca = 0.01 | 0.062 | 0.010 | 0.043 | 0.034
| Ca = 0.025 | 0.088 | 0.049 | 0.051 | 0.049
$d_{p}/d$ = 2/3 | Ca = 0.05 | 0.161 | 0.058 | 0.062 | 0.056
| Ca = 0.075 | Damaged | 0.220 | 0.105 | 0.101
| Ca = 0.1 | Damaged | 0.257 | 0.143 | 0.137
| Ca = 0.01 | 0.212 | 0.096 | 0.092 | 0.087
| Ca = 0.025 | 0.459 | 0.120 | 0.126 | 0.111
$d_{p}/d$ = 5/6 | Ca = 0.05 | Damaged | 0.372 | 0.159 | 0.158
| Ca = 0.075 | Damaged | Damaged | 0.228 | 0.224
| Ca = 0.1 | Damaged | Damaged | 0.317 | 0.302
Table 1: Measured surface strain $\delta S$ achieved by the capsule at the
impinging instant as a function of $d_{p}/d$ and Ca for $\theta$ = $\pi/6$
(both junction geometries), $\pi/3$ and $\pi/2$.

The deformation of the capsule navigating the bifurcation for different
occlusion rate is provided also for $\theta$ = $\pi/3$ and $\pi/2$. The
distributions of $\delta S$ over time as a function of $d_{p}/d$ and Ca are
collected in Supp.Figure.3. Interestingly, for larger $\theta$ almost
negligible deformations (smaller than 0.05) are registered for $d_{p}/d\leq
1/2$. For larger occlusion ratio, the trends enlighten for $\theta=\pi/6$ are
confirmed for large semi-aperture angles. For completeness, the measured
values of $\delta S$ at the impinging instant are documented in Table.1 for
all of the combinations of the investigated parameters.

###  2.3 Transport of multiple capsules in a bifurcation

A train of three initially spherical capsules with occlusion ratio $d_{p}/d$ =
2/3 is considered to further analyze the deformation of capsules crossing the
bifurcation. Specifically, two additional capsules are placed on top of the
first one, kept in the same position as before, namely 1.16 d from the inlet
section, with a mutual distance between the three equal to 0.85 d (see
Figure.10).

Figure 10: Train of three capsules navigating a microfluidic bifurcation. (a,
b, c) Trajectory of the capsules transported into the bifurcation for
$\theta=\pi/6$ (a), $\pi/3$ (b), and $\pi/2$ (c) and Ca = 0.01. (d, e, f)
Distribution of the relative stretching of the capsules for $\theta=\pi/6$
(d), $\pi/3$ (e), and $\pi/2$ (f) and Ca = 0.01. (h, g, i) Distribution of the
relative stretching of the capsules for $\theta=\pi/6$ (h), $\pi/3$ (g), and
$\pi/2$ (i) and Ca = 0.05. (l, m, n) Distribution of the relative stretching
of the capsules for $\theta=\pi/6$ (l), $\pi/3$ (m), and $\pi/2$ (n) and Ca =
0.1. Colors refer to the specific capsules.

Interestingly, when transporting multiple capsules into the bifurcation the
trajectory patterns complicate. When crossing the bifurcation all particles
may navigate one of the two branches in the same fashion as before. However,
as already observed, the total cross-section of the main channel is conserved
across the bifurcation so that the capsules have to squeeze to enter into a
daughter branch after crossing the junction region. Therefore, the established
flow results to be partially occluded by the presence of a capsule within one
of the two smaller channels, and the last particle (colored in red in
Figure.10) cannot be easily dragged across the junction zone. For Ca = 0.01
this phenomenon is documented in terms of traveled path in Figure.10a–.10c for
$\theta=\pi/6$, $\pi/3$, and $\pi/2$ while the entire capsules’ journey is
reported in Supp.movie.1–.3. Moreover, the distribution of the relative
stretching registered for the three capsules is compared with the one obtained
for the transport of a single capsule used as reference. Large stretching is
observed almost for all of the investigated combinations of $\theta$ and Ca
(see Figure.10d–.10n) while for Ca = 0.1 and $\theta=\pi/6$ the degradation of
both upper capsules is registered. Indeed, the relative stretching
distributions indicate that for $\theta=\pi/3$ and $\pi/2$ the second capsule
squeezes more than the other two for all Ca due to the presence of the first
capsule and the relative obstruction seen when crossing the junction zone (see
Supp.Movie.4-.6).

Figure 11: Train of rigid and soft capsules navigating the bifurcation for
$\theta=\pi/6$. (a, b, c) Trajectory of the capsules transported into the
bifurcation for $\theta=\pi/6$ when considering a rigid capsule in the top
(a), middle (b), and bottom (c) releasing position. d. Configuration of the
transported capsules for $\theta=\pi/6$ and with the rigid capsule in the
bottom releasing position from three different point of views for $tu_{max}/d$
= 0, 5, 15, 40.

The transport of microvesicles is assessed also when considering rigid objects
and their interaction with deformable membranes for $\theta=\pi/6$ and Ca =
0.01. In particular, three more configurations have been proposed by placing a
rigid plain sphere with the same density of soft capsules and diameter
$\frac{1}{2}d_{p}$ in place of a deformable capsule (a releasing position per
time). All of the initial conditions are reported in Figure.11a.

Figure 12: Train of three rigid capsules navigating the bifurcation for
$\theta=\pi/6$. a. Trajectory of the capsules transported into the
bifurcation. b. Capsules configurations during the transport from three
different point of views for $tu_{max}/d$ = 0, 5, 10, 40. c. Distribution of
the z component of the center of mass velocity for the three particle with a
close-up corresponding to the crossing of the junction region.

The transport of such objects along with their deformation and trajectories is
complicated again due to two phenomena: i) rigid objects would travel the
bifurcation faster than soft capsules since they are not subject to any
deformation crossing the junction region; ii) due to their size, rigid
particles would tend to avoid direct interaction with deformable membranes by
moving laterally and continuing to follow flow streamlines. The registered
trajectories of these three soft or rigid capsules are reported in
Figure.11a–.11c (see also Supp.Movie.7–.9). When the rigid particle is in the
upper releasing position it will travel across the bifurcation almost
unperturbed by the two following larger capsules (see Figure.11a). Note that,
again the last capsule (colored in red) follows the path chosen by the soft
particle above (colored in blue) due to the occurrence of the flow obstruction
and consequent capillary suction operated by the latter (see Supp.Movie.7). On
the contrary, when placing the rigid sphere in the bottom releasing position
(see Figure.11c), the latter need to flow laterally to avoid contact with the
two larger membranes in the main channel and then moves in a dislodged line
with respect to the duct axis trying into the selected daughter-branch. The
capsules configurations during these two moments are depicted in Figure.11d,
the small rigid sphere tends to marginate due to the interaction with larger
objects. This phenomenon is enhanced when transporting three rigid particles
within the microchip (see Figure.12). Such rigid particles do not interact
between them during the journey dislodging one with respect to the others to
find the smoother and more undisturbed path to overcome the bifurcation. For
completeness, also in this last case, the three particles’ journey is reported
in Supp.Movie.10. Moreover, the distribution of the streamwise component of
the center of mass velocity for the three particles is documented in
Figure.12c. The time interval in which the particle traverse the junction
region is enlightened into a close-up. The upper particle impinges against the
junction edge at 9.28 $tu_{max}/d$ (black line) and crawls on the right
branch. In the meantime, the immediately following particle (colored in blue)
selects the left branch thus avoiding any interaction with particles or with
the vasculature. For the same reason, the last particle (colored in red)
decides to follow the black particle that is already sufficiently far.

## Conclusions and future work

The deformation of inertial capsules navigating a microfluidic bifurcation as
a function of their stiffness and size is critically analyzed in this paper.
The incompressible Navier-Stokes equation is modeled through a lattice
Boltzmann scheme in which body forces account for immersed geometries. The
latter is treated with two different approaches; on one side, no-slip boundary
conditions are imposed on moving deformable or rigid geometries through the
dynamic-Immersed Boundary method, while, on the other side, on fixed immersed
geometries the Bouzidi-Firdaouss-Lallemand second-order bounce back technique
is adopted. As per the external boundaries, Dirichlet boundary conditions are
imposed with the Zou-He formulation.

The proposed model was adopted to systematically analyze the transport of
micrometric capsules when transported in a microfluidic bifurcation as a
function of their mechanical stiffness, size, bifurcation aperture angle,
main-to-daughter channels junction geometry, and the number of traveling
particles. In particular, initially spherical capsules with different
mechanical stiffness (Ca = 0.01, 0.025, 0.05, 0.075, and 0.1) and different
size ($d_{p}/d=$ 1/3, 1/2, 2/3, and 5/6) are transported one per time within a
bifurcation with semi-aperture $\theta=\pi/6$. This bifurcation consists of
the main channel that splits into two daughter ducts. Such splitting is
performed with two junctions: a sharp joining of the three channels and a more
smooth joint with a trapezoidal vertical section. For $d_{p}/d<$ 2/3 the
smooth junction is found to not injure the membranes transported within thanks
to the flat edge on top of the junction connecting the two branches. On the
contrary, the sharp junction is found to imply strong deformations on the
navigating capsules. Large particle ($d_{p}/d$ = 5/6) are damaged with both
junctions for Ca $>$ 0.05. Moreover, it is demonstrated that for larger
aperture angles ($\theta=$ $\pi/3$ and $\pi/2$) smaller deformations are
registered for all particles size. The measured distribution of the relative
stretching during the transport is documented for all apertures as a function
of Ca and $d_{p}/d$ demonstrating that small particles ($d_{p}/d$ = 1/3) exert
almost null deformations while larger capsules ($d_{p}/d\geq$ 2/3) may be
damaged also with the smoother junction for large values of Ca. Then, the
transport of a train of three capsules into the bifurcation is assessed for
$\theta=\pi/6$ as a function of Ca. Interestingly, when multiple capsules
navigate the bifurcation the particle released into the lowest releasing
position deal with the flow obstruction due to the presence of the other two
membranes. After one of the two obstructing capsules reach the bifurcation
peripheries, the flow is restored and the last particle is eventually dragged
into a daughter branch. Lastly, those three deformable capsules are
substituted with a rigid smaller particle one per time. For the specific set
of investigated parameters, smaller rigid particles are found able to dislodge
laterally in order to avoid direct interactions with larger capsules.

Collectively, these data provide specific insights for the scientific
community investing efforts in the rational design of microfluidic chips for
the assessment of particle dynamics and cell sorting in human
microcirculation. The spanned region of the phase-space determined by the
number of the investigated parameters deepens the knowledge of the fundamental
physics behind the prototypical problem of a particle navigating a
bifurcation.

## Acknowledgments

AC acknowledges the project “Research for Innovation” (REFIN) - POR Puglia
FESR FSE 2014-2020 - Asse X - Azione 10.4 (Grant No. CUP - D94120001410008).
AC and TP are members of Gruppo Nazionale per il Calcolo Scientifico (GNCS) of
Istituto Nazionale di Alta Matematica (INdAM). This research has been
partially supported by the Research Project of National Relevance “Multiscale
Innovative Materials and Structures” granted by the Italian Ministry of
Education, University and Research (MIUR Prin 2017, project code 2017J4EAYB
and the Italian Ministry of Education, University and Research under the
Programme Department of Excellence Legge 232/2016 (Grant No. CUP -
D94I18000260001).

## Competing Interests

The authors declare no competing interests.

## References

  * Guidolin et al. [2011] Diego Guidolin, Enrico Crivellato, and Domenico Ribatti.  The “self-similarity logic” applied to the development of the vascular system.  _Developmental biology_ , 351(1):156–162, 2011\.
  * Stoecklein and Di Carlo [2019] Daniel Stoecklein and Dino Di Carlo.  Nonlinear microfluidics.  _Anal. Chem._ , 91(1):296–314, 2019.
  * Volpe et al. [2017] A Volpe, P Paiè, A Ancona, R Osellame, P M Lugarà, and G Pascazio.  A computational approach to the characterization of a microfluidic device for continuous size-based inertial sorting.  _Journal of Physics D: Applied Physics_ , 50(25):255601, 2017.
  * Koolivand and Dimitrakopoulos [2017] A Koolivand and P Dimitrakopoulos.  Deformation of an elastic capsule in a microfluidic t-junction: settling shape and moduli determination.  _Microfluidics and Nanofluidics_ , 21(5):89, 2017\.
  * Gandhi et al. [2020] Tanvi Gandhi, Jinzi Mac Huang, Antoine Aubret, Yaocheng Li, Sophie Ramananarivo, Massimo Vergassola, and Jérémie Palacci.  Decision-making at a t-junction by gradient-sensing microscopic agents.  _Physical Review Fluids_ , 5(10):104202, 2020\.
  * Tran et al. [2020] SBQ Tran, QT Le, FY Leong, and DV Le.  Modeling deformable capsules in viscous flow using immersed boundary method.  _Physics of Fluids_ , 32(9):093602, 2020.
  * Doyeux et al. [2011] Vincent Doyeux, Thomas Podgorski, Sarah Peponas, Mourad Ismail, and Gwennou Coupier.  Spheres in the vicinity of a bifurcation: elucidating the zweifach–fung effect.  _Journal of Fluid Mechanics_ , 674:359–388, 2011.
  * Sevenie et al. [2015] Benjamin Sevenie, Vincent Fremont, Dominique Barthès-Biesel, and Anne-Virginie Salsac.  Flow of an artificial capsule in a bifurcation: from image processing to shape modal description.  In _Congrès français de mécanique_. AFM, Association Française de Mécanique, 2015.
  * Wang et al. [2018] Z Wang, Y Sui, A-V Salsac, D Barthès-Biesel, and W Wang.  Path selection of a spherical capsule in a microfluidic branched channel: Towards the design of an enrichment device.  _Journal of Fluid Mechanics_ , 849:136–162, 2018.
  * Coclite and Gambaruto [2019] Alessandro Coclite and Alberto M. Gambaruto.  Injection of deformable capsules in a reservoir: A systematic analysis.  _Fluids_ , 4(3), 2019.  ISSN 2311-5521.  doi: 10.3390/fluids4030122.
  * Lu et al. [2021] R. X. Lu, Z. Wang, A.-V. Salsac, D Barthes-Biesel, W Wang, and Y Sui.  Path selection of a train of spherical capsules in a branched microchannel.  _Journal of Fluid Mechanics_ , 923:A11, 2021.
  * Hymel et al. [2019] Scott J Hymel, Hongzhi Lan, Hideki Fujioka, and Damir B Khismatullin.  Cell trapping in y-junction microchannels: A numerical study of the bifurcation angle effect in inertial microfluidics.  _Physics of Fluids_ , 31(8):082003, 2019.
  * Häner et al. [2020] Edgar Häner, Matthias Heil, and Anne Juel.  Deformation and sorting of capsules in a t-junction.  _Journal of Fluid Mechanics_ , 885, 2020.
  * Shields IV et al. [2015] C Wyatt Shields IV, Catherine D Reyes, and Gabriel P López.  Microfluidic cell sorting: a review of the advances in the separation of cells from debulking to rare cell isolation.  _Lab on a Chip_ , 15(5):1230–1249, 2015.
  * Vesperini et al. [2017] Doriane Vesperini, Oriane Chaput, Nadège Munier, Pauline Maire, Florence Edwards-Levy, Anne-Virginie Salsac, and Anne Le Goff.  Deformability-and size-based microcapsule sorting.  _Medical Engineering & Physics_, 48:68–74, 2017.
  * Häner et al. [2021] Edgar Häner, Doriane Vesperini, Anne-Virginie Salsac, Anne Le Goff, and Anne Juel.  Sorting of capsules according to their stiffness: from principle to application.  _Soft Matter_ , 17(13):3722–3732, 2021.
  * Guzniczak et al. [2020] Ewa Guzniczak, Oliver Otto, Graeme Whyte, Nicholas Willoughby, Melanie Jimenez, and Helen Bridle.  Deformability-induced lift force in spiral microchannels for cell separation.  _Lab on a Chip_ , 20(3):614–625, 2020.
  * Manneschi et al. [2016] C Manneschi, RC Pereira, G Marinaro, A Bosca, M Francardi, and P Decuzzi.  A microfluidic platform with permeable walls for the analysis of vascular and extravascular mass transport.  _Microfluidics and Nanofluidics_ , 20(8):1–12, 2016.
  * Miali et al. [2019] Marco E Miali, Marianna Colasuonno, Salvatore Surdo, Roberto Palomba, Rui Pereira, Eliana Rondanina, Alberto Diaspro, Giuseppe Pascazio, and Paolo Decuzzi.  Leaf-inspired authentically complex microvascular networks for deciphering biological transport process.  _ACS applied materials & interfaces_, 11(35):31627–31637, 2019.
  * Mollica et al. [2019] Hilaria Mollica, Roberto Palomba, Rosita Primavera, and Paolo Decuzzi.  Two-channel compartmentalized microfluidic chip for real-time monitoring of the metastatic cascade.  _ACS Biomaterials Science & Engineering_, 5(9):4834–4843, 2019.
  * Melchionna [2011] Simone Melchionna.  A model for red blood cells in simulations of large-scale blood flows.  _Macromolecular Theory and Simulations_ , 20(7):548–561, 2011.
  * Freund [2014] Jonathan B. Freund.  Numerical simulation of flowing blood cells.  _Ann. Rev. Fliud Mech._ , 46:67–95, 2014.
  * Ye et al. [2014] Swe Soe Ye, Yan Cheng Ng, Justin Tan, Hwa Liang Leo, and Sangho Kim.  Two-dimensional strain-hardening membrane model for large deformation behavior of multiple red blood cells in high shear conditions.  _Theoretical Biology and Medical Modelling_ , 11(1):19, 2014.
  * Secomb [2017] Timothy W. Secomb.  Blood flow in the microcirculation.  _Ann. Rev. Fliud Mech._ , 49:443–461, 2017.
  * Sebastian and Dittrich [2018] Bernhard Sebastian and Petra S. Dittrich.  Microfluidics to mimic blood flow in health and disease.  _Ann. Rev. Fliud Mech._ , 50:483–504, 2018.
  * Laurino et al. [2020] Federica Laurino, Alessandro Coclite, Annagiulia Tiozzo, Paolo Decuzzi, and Paolo Zunino.  A hierarchical multiscale model for predicting the vascular behavior of blood-borne nanomedicines.  _International Journal for Multiscale Computational Engineering_ , 18(3), 2020.  doi: 10.1615/IntJMultCompEng.2020033358.
  * Coclite [2020] Alessandro Coclite.  Vascular journey and adhesion mechanics of micro-sized carriers in narrow capillaries.  _Microvascular Research_ , 132:104069, 2020.  ISSN 0026-2862.  doi: https://doi.org/10.1016/j.mvr.2020.104069.
  * Coclite et al. [2019] A. Coclite, S. Ranaldo, M.D. de Tullio, P. Decuzzi, and G. Pascazio.  Kinematic and dynamic forcing strategies for predicting the transport of inertial capsules via a combined lattice Boltzmann immersed boundary method.  _Computers & Fluids_, 180:41–53, 2019.  ISSN 0045-7930.  doi: https://doi.org/10.1016/j.compfluid.2018.12.014.
  * Coclite et al. [2020] Alessandro Coclite, Sergio Ranaldo, Giuseppe Pascazio, and Marco D. de Tullio.  A lattice boltzmann dynamic-immersed boundary scheme for the transport of deformable inertial capsules in low-re flows.  _Computers & Mathematics with Applications_, 80(12):2860 – 2876, 2020.  ISSN 0898-1221.  doi: https://doi.org/10.1016/j.camwa.2020.09.017.
  * Zou and He [1997] Qisu Zou and Xiaoyi He.  On pressure and velocity boundary conditions for the lattice Boltzmann BGK model.  _Physics of Fluids_ , 9(6):1591–1598, 1997.  doi: http://dx.doi.org/10.1063/1.869307.
  * Bouzidi et al. [2001] M’hamed Bouzidi, Mouaouia Firdaouss, and Pierre Lallemand.  Momentum transfer of a boltzmann-lattice fluid with boundaries.  _Physics of fluids_ , 13(11):3452–3459, 2001\.
  * Bhatnagar et al. [1954] P. L. Bhatnagar, E. P. Gross, and M. Krook.  A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems.  _Phys. Rev._ , 94:511–525, May 1954.  doi: 10.1103/PhysRev.94.511.
  * Guo et al. [2002] Zhaoli Guo, Chuguang Zheng, and Baochang Shi.  Discrete lattice effects on the forcing term in the lattice Boltzmann method.  _Phys. Rev. E_ , 65:046308, Apr 2002.  doi: 10.1103/PhysRevE.65.046308.
  * Coclite et al. [2014] A. Coclite, G. Gonnella, and A. Lamura.  Pattern formation in liquid-vapor systems under periodic potential and shear.  _Phys. Rev. E_ , 89:063303, Jun 2014.  doi: 10.1103/PhysRevE.89.063303.
  * Coclite et al. [2016] A. Coclite, M. D. de Tullio, G. Pascazio, and P. Decuzzi.  A combined lattice Boltzmann and immersed boundary approach for predicting the vascular transport of differently shaped particles.  _Computers & Fluids_, 136:260 – 271, 2016.  ISSN 0045-7930.  doi: http://dx.doi.org/10.1016/j.compfluid.2016.06.014.
  * Peskin [2002] Charles S. Peskin.  The immersed boundary method.  _Acta Numerica_ , 11:479–517, 1 2002.  doi: 10.1017/S0962492902000077.
  * de Tullio and Pascazio [2016] M. D. de Tullio and G. Pascazio.  A moving-least-squares immersed boundary method for simulating the fluid-structure interaction of elastic bodies with arbitrary thickness.  _Journal of Computational Physics_ , submitted, 2016.
  * Vanella and Balaras [2009] Marcos Vanella and Elias Balaras.  A moving-least-squares reconstruction for embedded-boundary formulations.  _Journal of Computational Physics_ , 228(18):6617 – 6628, 2009.  ISSN 0021-9991.  doi: http://dx.doi.org/10.1016/j.jcp.2009.06.003.
  * Gelder [1998] Allen Van Gelder.  Approximate simulation of elastic membranes by triangulated spring meshes.  _Journal of graphics tools_ , 3(2):21–41, 1998\.
  * Tan et al. [2012] J. Tan, A. Thomas, and Y. Liu.  Influence of red blood cells on nanoparticle targeted delivery in microcirculation.  _Soft Matter_ , 8(6):1934–1946, 2012.  doi: 10.1039/c2sm06391c.