Dissipative Particle Dynamics Study of Strain Distribution in Capsules Deformed by Microfluidic Constrictions

Nishanthi Rajkamal¹ and Srikanth Vedantam^2,*

¹Department of Biotechnology, Indian Institute of Technology Madras, Chennai 600036, India
²Department of Engineering Design, Indian Institute of Technology Madras, Chennai 600036, India
E-mail: nishanthir169@gmail.com; srikanth@iitm.ac.in
${}^{*}$Corresponding Author

Received 27 July 2021; Accepted 04 November 2021; Publication 27 November 2021

Abstract

We present a dissipative particle dynamics (DPD) study of the deformation of capsules in microchannels. The strain in the membrane during this deformation causes the formation of temporary pores, which is termed mechanoporation. Mechanoporation is being considered as a means by which intracellular delivery of a broad range of cargo can be facilitated. In this work, we examine the strain distribution on the capsule membrane during transport of the capsule in converging-diverging microchannels of different constriction widths. The pore density is correlated to the strain in the membrane. We find that the highest strains and, consequently, the highest pore densities occur at intermediate channel widths. This occurs due to a competition of the bending of the membrane and fluid shear stresses in the flow.

Keywords: Intracellular delivery, mechanoporation, dissipative particle dynamics, numerical simulation, microfluidics modeling.

1 Introduction

Intracellular delivery of macromolecules is an important step in research and therapeutics. For this, the cell membrane is disrupted to create pores or pore-like openings. Membrane disruption can be done by microinjection, electroporation, optoporation, a thermally induced process, biochemical agents or mechanoporation [1]. Fluid shear stresses in microfluidic channels exert mechanical forces on the cell membrane to cause mechanoporation. Industrial applications of microfluidic devices using constrictions in microfabricated channels are being developed for many cell-based therapies [2–8]. For example, Sharei et al. [7] present a study of cell squeezing in a parallel network of microfluidic channels for intracellular delivery. The cells experience strain in the membrane due to the deformation in the channels. The strain in the cell membrane during transport in the microchannels disrupts the plasma membrane to create pores which facilitate intracellular delivery. The pores in the membrane are mechanosensitive Stretch Activated (SA) channels, which open when the cell membrane is stretched by mechanical forces. Macromolecules in the suspension diffuse through these SA channels into the cell. The plasma membrane undergoes recovery by closing these openings upon unloading [9].

Higher strain on the cell during deformation may cause off-target damages to the intracellular structures such as the cytoskeleton, nucleus and the genomic DNA [10]. These conditions stimulate cell death pathways [11] which, in turn, lead to low cell viability and reduce device efficiency. Experiments [7] have shown that the channel geometry parameters such as the length, width and configuration of the constrictions affect the cell viability and device efficiency. Hence a detailed study of the strain experienced by cells during deformation in the constriction is essential for developing devices with minimal cell damage during the process.

Pak et al. [12] recently developed a multiscale continuum approach to model gating of SA channels. They considered 2D vesicles being transported in constriction channels. The dynamic opening and closing of the SA channels was modeled based on membrane tension by placing the channels in one of 64 locations. The channel was considered to be open only when the membrane tension is sufficiently large, otherwise it remained closed. The channels in open state were only found to be in the front end of the vesicle. Since the SA channel did not move around the vesicle membrane, the initial position was important in determination of the open or closed states [12].

However, it is important to note that the strain distribution in the vesicle membrane is highly non-uniform and is affected by the microfluidic channel geometry and fluid flow distribution. In order to model these effects in greater detail we propose an alternative approach. We model the fluid flow using dissipative particle dynamics (DPD) [13] and the capsule by means of discrete particles interacting through various interparticle forces. Through this, we are able to model the distribution of strain around the membrane and correlate it with the pore density. This approach is more general than studying the effect on channels in specified locations.

The study of hydrodynamics of fluids in the macroscopic scale is done by considering the fluid as a continuum domain using computational fluid dynamics (CFD). When the fluids are examined at smaller length scales the hydrodynamics observed at the macroscopic scale still holds good but the discreteness of the fluid manifests itself in the form of Brownian effects. These Brownian effects and other van der Waals forces also play a significant effect on the hydrodynamics and interactions with particles. This length scale between the macroscopic region and the atomistic scale called as the mesoscale region is in the range of 10–1000 $nm$ and 1 $ns$ – 10 $\mu s$. Molecular dynamics (MD) simulations which deals with the interaction between atoms is one way to study these interactions but it is computationally expensive. Obtaining hydrodynamics at the mesocale through MD simulations would require extremely large number of molecules $(\sim {10}^{21})$ and time scale much larger than collision time $({10}^{-15}s)$. Over the last three decades several mesoscopic techniques have been developed to deal with the hydrodynamics by neglecting some of the excessive short time and length scale details as seen in MD. Some of the widely researched methods are the Lattice-Boltzmann method (LBM) [14], Brownian dynamics (BD), and lattice-gas cellular automata (LGA) [15].

There have been several mesoscopic models developed to study the deformable nature of RBCs, like the Dzwinel et al. model of RBCs as an elastic material [16], the Lattice Boltzmann method [17] and the Multiparticle Collision Dynamics [18]. Pivkin et al. used DPD to create an accurate coarse-grained model to study RBCs with single elastic spring networks [19]. This was further developed by Fedosov et al. [20] by developing a network of viscoelastic springs. Later, a low-dimensional model was developed by Pan et al. which was computationally efficient and also captured RBC dynamics in vessels with large diameter [21]. A detailed comparison and qualitative analysis of the 2D and 3D model for the blood components has been studied in detail by Müller [22]. The similarities in margination, haematocrit and shear rate of RBCs, platelets and proteins in the high and low dimension models were established. Deformation and dynamics of RBCs in shear flows and confinement flows have also been studied in both 2D [23] and 3D [20]. The difference in deformation of healthy and infected RBC when flowing in stenosed channel based on the cell stiffness have been examined [23]. There have been several studies using similar 2D models to study the dynamics of an elastic capsules, vesicles and cells [24–27].

The shape transitions of the initial spherical capsules into bullet-like shapes and parachute-like shapes based on the confinement of flow and the shear forces have been studied widely in experiments and computational simulations in various channel geometries [28–31]. Several shapes observed in experiments [30] were captured by the 2D capsule models [32]. In this work, our focus of interest is the strain distribution on the membrane due to different deformation shapes of the cells. These have been shown to be qualitatively captured by the 2D models [32].

In this paper, we present a DPD approach to simulate hydrodynamics of capsule transported in constriction channels. We study changes in strain experienced in the membrane for different constriction widths. The strain on the cell due to fluid shear as well as deformation is studied using this model. Competition between the bending of membrane due to deformation and the fluid shear stress affects the strain distribution. We use experimental data [33] to correlate strain with pore density. We find that the channels of intermediate width has the highest strain and consequently, highest pore density.

2 DPD Model

2.1 Model of the Fluid

In this section, we briefly outline the DPD model [13], particularly, the Finite sized DPD (FDPD) model proposed by Pan et al. [34]. The fluid is taken to consist of $N$ uniformly sized particles of radius $R$. The number density of the particles per unit area is denoted by $\rho$. Since the FDPD particles are taken to be of uniform size, their masses and mass moments of inertia are all the same and denoted by $m$ and $I$ respectively. The particles obey the laws of conservation of linear and angular momentum,

$m{\displaystyle \frac{\mathrm{d}{\mathbf{\text{v}}}_i}{dt}}$	$={\displaystyle \sum _{j\ne i}}{\mathbf{\text{f}}}_{ij},$	(1)
$I{\displaystyle \frac{\mathrm{d}{\mathit{\omega }}_i}{dt}}$	$={\displaystyle \frac{1}{2}}{\displaystyle \sum _{j\ne i}}{\mathit{r}}_{ij}\times {\mathbf{\text{f}}}_{ij}$	(2)

where ${\mathbf{v}}_i$ and ${\mathit{\omega }}_i$ are the velocity and angular velocity of the $i^{th}$ particle respectively and ${\mathit{r}}_{ij}={\mathit{r}}_i-{\mathit{r}}_j$. The forces ${\mathbf{f}}_i$ arise from the interaction of the $i^{th}$ particle with all of its neighbors $j$. Thus the sums in (1) and (2) are taken over all particles $j$. Practically, the interactions are limited to a neighborhood of radius $r_c$ around the $i^{th}$ particle. The torque due to the interaction of particle $j$ on particle $i$ is taken around the mid-point of the line joining the centers of the two particles.

The interaction forces play a key role in the behavior of this model. As is standard, the interaction forces ${\mathbf{f}}_{ij}$ comprise of four kinds of forces: (i) conservative $({\mathbf{f}}^C)$, (ii) dissipative central $({\mathbf{f}}^T)$, (iii) dissipative rotational $({\mathbf{f}}^R)$ and (iv) stochastic $({\mathbf{f}}^S)$. Thus the total force on the particle $i$ due to its neighbor $j$ is given by,

$${\mathbf{f}}_{ij}={\mathbf{f}}_{ij}^C+{\mathbf{f}}_{ij}^T+{\mathbf{f}}_{ij}^R+{\mathbf{f}}_{ij}^S$$

(3)

and the total force on $i^{th}$ particle due to all its neighbours is given by,

$${\mathbf{f}}_i=\sum _{j\ne i}{\mathbf{f}}_{ij}+{\mathbf{f}}^{\mathrm{E}}$$

(4)

The conservative force is a soft sphere interaction acting along the line connecting the particle centers given by

$${\mathbf{f}}_{ij}^{\mathrm{C}}=a_{ij}\mathrm{\Gamma }(r_{ij}){\widehat{\mathbf{\text{e}}}}_{ij},$$

(5)

where $a_{ij}$ is stiffness, $r_{ij}=|{\mathit{r}}_{ij}|$ and ${\widehat{\mathbf{\text{e}}}}_{ij}={\mathit{r}}_{ij}/r_{ij}$ is a unit vector. An appropriate weight function, $\mathrm{\Gamma }(r_{ij})$ is selected so that the conservative force vanishes at the cutoff radius $r_{ij}=r_c$. It is common to take

$$\mathrm{\Gamma }(r_{ij})=\{\begin{array}{cc}1-\frac{r_{ij}}{r_c},\hfill & \text{for}{r}_{ij}<r_c\hfill \\ 0\mathit{\hspace{1em}}\hfill & \text{for}{r}_{ij}>r_c.\hfill \end{array}$$

(6)

The translation force is assumed to have central and non-central dissipative components given by

${\mathbf{f}}_{ij}^T=-{\gamma }_{ij}^{||}{\mathrm{\Gamma }}^2(r_{ij})({\mathbf{v}}_{ij}\cdot {\widehat{\mathbf{e}}}_{ij}){\widehat{\mathbf{\text{e}}}}_{ij}-{\gamma }_{ij}^{\perp }{\mathrm{\Gamma }}^2(r_{ij})({\mathbf{v}}_{ij}-({\mathbf{v}}_{ij}\cdot {\widehat{\mathbf{\text{e}}}}_{ij}){\widehat{\mathbf{\text{e}}}}_{ij})$

(7)

where ${\gamma }_{ij}^{||}$ and ${\gamma }_{ij}^{\perp }$ are longitudinal and transverse dissipation coefficients. This dissipative force reduces the relative velocity ${\mathbf{v}}_{ij}={\mathbf{v}}_i-{\mathbf{v}}_j$ between particles in both directions. The rotational dissipative force is taken to be of the form

$${\mathbf{f}}_{ij}^R=-{\gamma }_{ij}^{\perp }{\mathrm{\Gamma }}^2(r_{ij})\left({\mathit{r}}_{ij}\times \frac{1}{2}({\mathit{\omega }}_i+{\mathit{\omega }}_j)\right)$$

(8)

A stochastic force accounting for the lost degrees of freedom of individual molecules due to coarse-graining is taken to be

${\mathbf{f}}_{ij}^S\mathrm{\Delta }t=\mathrm{\Gamma }(r_{ij})\left\{{\sigma }_{ij}^{||}tr[d{\mathbf{W}}_{ij}]{\displaystyle \frac{1}{\sqrt{d}}}\mathbf{1}+\sqrt{2}{\sigma }_{ij}^{\perp }d{\mathbf{W}}_{ij}^A\right\}\cdot {\widehat{\mathbf{\text{e}}}}_{ij}$

(9)

where $\mathrm{\Delta }t$ is the time step, $d=2$ for two-dimensional simulations and $tr[d{\mathbf{W}}_{ij}]$ is trace of the symmetric independent Wiener increment matrix $d{\mathbf{W}}_{ij}$, while $d{\mathbf{W}}_{ij}^A$ is its antisymmetric part. According to the fluctuation dissipation theorem, the random and dissipation coefficients are correlated by ${\sigma }_{ij}^{\parallel }=\sqrt{2k_{B}T{\gamma }_{ij}^{\parallel }}$ and ${\sigma }_{ij}^{\perp }=\sqrt{2k_{B}T{\gamma }_{ij}^{\perp }}$. The stochastic and dissipation forces together maintain a constant granular temperature of the system.

2.2 Model of the Capsule

The capsule is modelled as a two dimensional closed bead-spring chain [22, 35] as shown in Figure 1. The particles interact through elastic interactions with the potential energy given by

$$E_{\mathit{Spring}}=k_{B}T\frac{k_S}{2}\sum _{\mathrm{i}=1}^{N_c}{\left(\frac{l_i-l_0}{l_0}\right)}^2$$

(10)

where $l_i$ is the length of the $i^{th}$ spring at a given instant and $l_0$ is the reference length of each of the springs. The spring constant $k_S$ represents the resistance to membrane stretching deformation of the capsule. In order to maintain a circular shape radial springs with spring constant $k_R$ is considered to connect the particle on the membrane to a particle in the centre. Finally, an area conservation constraint is also considered

$$E_{\mathit{Area}}=k_{B}T\frac{k_A}{2}\sum _{\mathrm{i}=1}^{N_c}{\left(\frac{A_t-A_0}{A_0}\right)}^2$$

(11)

where $k_A$ is the area constraint constant, $A_t$ is the area of the capsule at any instant and $A_0$ is the initial total area. Generally, $A_0$ is taken to be the area of a circle with the desired diameter of the capsule. The diameter of the capsule $D_c$ is chosen to be 12 $\mu m$ corresponding to that of a HeLa cell [7].

Figure 1 Representation of a capsule. The capsule is taken to be composed of several DPD particles on the boundary and a central particle interacting with them along the radial lines.

The interaction forces due to neighboring particles on $i^{th}$ particle is given by

$${\mathbf{f}}_i^S=-\frac{\partial E_{\mathit{Spring}}}{\partial {\mathbf{r}}_i},{\mathbf{f}}_i^A=-\frac{\partial E_{\mathit{Area}}}{\partial {\mathbf{r}}_i}$$

(12)

2.3 Boundary Conditions

In order to impose a no-slip boundary condition on the channel surfaces we employ the approach outlined in Ranjith et al. [36]. When a fluid particle is within a cut-off distance $(r_c)$ from the wall, it is taken to interact with a DPD particle on the nearest point at the wall. In other words the DPD particle with the coordinates $(x_p,y_p)$ at a distance less than $r_c$ interacts with the wall at the closest point $(x_w,y_w)$ which is given by,

$x_w$	$={\displaystyle \frac{1}{(1+m^2)}}(x_p+m{y}_p+m{y}_b-m^2{x}_b)$	(13)
$y_w$	$={\displaystyle \frac{1}{(1+m^2)}}(m{x}_p+y_b-m{x}_b+m^2{y}_p)$	(14)

where $m$ is the slope of the inclined wall and $(x_b,y_b)$ is an end point of the inclined wall (Figure 2). This model of the wall boundary eliminates the crystalline structure of fixed particle approaches and has been shown to reduce spurious density fluctuations [36].

Figure 2 Schematic of boundary condition implementation. An IFP particle is identified on the wall at locations nearest to the DPD particles.

The no-slip condition is achieved by controlling the lateral dissipative force component between the fluid particle and the wall particle along the direction tangential to the wall. The lateral dissipation co-efficient is taken to be ${\gamma }_{\mathit{fw}}^{\perp }=\alpha {(1-r_{\mathit{fw}}/r_c)}^2{\gamma }_{\mathit{ff}}^{\perp }$ within the cut-off distance $r_c$ from the wall. A slip modification parameter $\alpha$ is used to model no-slip at the wall. Similar to fluid-fluid interactions the dissipative and the random coefficients are related through ${\sigma }_{\mathit{fw}}^{\perp }=\sqrt{2k_{B}T{\gamma }_{\mathit{fw}}^{\perp }}$. Periodic boundary conditions are applied at inlet and outlet region which implies a particle exiting at the outlet is re-introduced at the inlet with the same velocity. The capsule and wall interactions are modeled similar to fluid-wall interactions.

2.4 Model Parameters

The size of the two-dimensional simulation domain is taken to be $300r_c\times 40r_c$ with a number density of $\rho =3$. The x-axis is along the channel length which is also the direction of fluid flow and y-axis is along the channel width. The channel is divided into five regions namely – inlet region, converging region, constriction region, diverging region and outlet region as shown in Figure 3. The upstream region (inlet and converging regions) is prior to the constriction channel whereas the downstream region (outlet and diverging region) is the region following the constriction channel in the direction of fluid flow. In the experimental study of Sharei et al. [7], the focus of the channel design was mainly on the constriction region, the length of the upstream and downstream region and the slope of the converging and diverging region was not given much importance. The parameters here are chosen appropriately to avoid fluid recirculation zones and focusing flows. The length of the constriction region is $20r_c$. The other regions have an equal length of $70r_c$. The constriction region width, $W=12r_c$ is chosen, as it is approximately $50\%$ smaller than the diameter of the capsule. For analysis of strain on the membrane for varying constriction width, a narrower channel $W=8r_c$ and a wider channel $W=16r_c$ are studied.

Figure 3 Schematic of the channel showing different sections.

The capsule was constructed with $N_c=80$ particles and diameter $D_c=24r_c$ [35]. The parameters $k_S$ and $k_R$ are taken to be 800 and 1600 respectively while $k_A$ is 6000. These values ensure that the contour length and the area of the capsule would be constant while the capsule is being squeezed in the constriction region. The cut-off distance, $r_c=1.5$ [21], is for interaction within the capsule particles, whereas for rest of the interactions such as fluid-fluid, fluid-capsule, fluid-wall, capsule-wall, $r_c=1.0$. The temperature of the system is maintained constant at $k_{B}T=0.1$. The conservative force parameter is taken to be $a=50$ which is not inline with the relation $a=\frac{75k_{B}T}{\rho }$. The higher $a$ ensures the absence of pseudo-compressibility effects and a better particle distribution between upstream and downstream regions [37]. The same $a$ is used for conservative force interaction within the capsule particles. The dissipative force parameters ${\gamma }_{\mathit{ff}}^{\parallel }={\gamma }_{\mathit{ff}}^{\perp }={\gamma }_{\mathit{fw}}^{\parallel }=45.0$ and the coefficients of the stochastic force is ${\sigma }_{\mathit{ff}}^{\parallel }={\sigma }_{\mathit{ff}}^{\perp }={\sigma }_{\mathit{fw}}^{\parallel }=3.0$ [36]. The coefficient of the lateral dissipative component of fluid-wall and capsule-wall interactions differ as explained in Section 2.3.

A steady external force ${\mathbf{f}}^{\mathrm{E}}=0.06$ is applied on every particle along x-axis in the direction of flow. The temporal evolution of particle positions is obtained by integrating the momentum Equations (1) and (2) by implementing the modified velocity-Verlet computational scheme with a time step of $\mathrm{\Delta }t=0.005$.

3 Results and Discussion

3.1 Flow Through a Constriction Channel

The time-averaged velocity vector plot and the density contour plot are given in Figures 4(a) and 4(b) respectively. As explained earlier the higher $a$ and the chosen channel geometry aide in better particle distribution throughout the channel and absence of recirculation zones and focusing flows. The $Re$ numbers for these simulations are estimated to be in the range of $0.1$ to $1.0$.

Figure 4 A tapered microfluidic channel showing (a) time-averaged velocity vector plot and (b) density contour plot.

3.2 Capsule Mechanics

Different shapes of the deformed capsule are observed in plane Poiseuille flow at varying flow rates. At low flow rates the capsule appears to be circular and then deforms to form an egg shape as shown in Figures 5(a) and 5(b). As the flow rate is increased further the capsule forms a triangle shape Figure 5(c) and then an arrowhead or a bullet shape in Figure 5(d). The shapes are comparable to experimental findings of Hwang et al [30] as seen in Figures 5(e), (f), (g) and (h) for different shear rates.

Figure 5 Deformation of capsule for varying flow rate in plane poiseuille flow from simulations in the top and range of shapes observed for elastic particles from experiments by Hwang et al. [30] in the bottom. (a) circular, ${\mathbf{f}}^{\mathrm{E}}=0.04$, (b) egg, ${\mathbf{f}}^{\mathrm{E}}=0.08$, (c) triangular, ${\mathbf{f}}^{\mathrm{E}}=0.2$, (d) arrowhead or bullet, ${\mathbf{f}}^{\mathrm{E}}=0.5$, (e) circular, Capillary number, $Ca\approx 0.005$, (f) egg, $Ca\approx 0.05$, (g) triangular, $Ca\approx 0.1$, (h) arrowhead, $Ca\approx 0.2$. The scale bars in the bottom indicate 200 $\mu m$.

3.3 Strain on Capsule Flowing in Channel with Varying Constriction

As the capsule is transported in a channel of uniform width it assumes a classic egg shape as seen in Figure 6. The capsule is drawn at equal time intervals and it is observed that it is transported at constant speed. We note that in the steady state, there is an exchange of DPD particles across the capsule membrane but this does not affect the overall deformation characteristics.

Figure 6 Snapshots of a capsule transported in a poiseuille flow field at equal time intervals. The capsule deforms into a egg shape at steady state.

When the capsule is transported in a tapered channel, we observe an elongation of the capsule in the tapered section. As it exits the tapered section, it mushrooms due to the hydrodynamics of the fluid flow. The speed of the capsule also varies, with a lower velocity in the inlet and converging region of the channel and higher velocity at the exit of the constriction as seen in Figure 7. This is in accordance with hydrodynamics of flow in a tapered channel.

Figure 7 Snapshots of a capsule transported in a tapered channel at equal time intervals. The snapshots are close together at the entry of the narrow section indicating low speed whereas the speed increases at the exit of the constriction.

The strain $\epsilon$ in the capsule membrane is calculated as

$$\epsilon =\frac{(l-l_0)}{l_0}$$

(15)

where $l_0$ is the average length of the particles comprising the capsule membrane in a quiescent condition and $l$ is the time-averaged distance between the particles comprising the capsule membrane during transport in a constriction channel. A positive strain indicates that the membrane is under extension whereas a negative strain occurs where the membrane is under compression.

Figure 8 Deformation of the capsule membrane at the entry to the constriction in microchannels of varying constriction widths (a) W $=$ 8rc, (b) W $=$ 12rc and (c) W $=$ 16rc.

When the capsule is squeezed into a narrow channel, it develops a positive strain in some regions and negative strain in others as seen in Figure 8. The differences in the strain observed are due to the combined effect of shear stress and the deformation by the constriction. Though the extent of deformation of capsule is high in narrowest channel, $W=8r_c$, (Figure 8(a)), the strain observed is accommodated at the front of the capsule through its curvature. In the case of deformation of capsule in the intermediate width channel, $W=12r_c$, the curvature at the front is lower but the strain is found to be greater, owing to the additional strain caused by shear effects at the wall as seen in Figure 8(b). In the widest channel, $W=16r_c$, (Figure 8(c)), the deformation is the least due to lower shear stresses of the wall.

Figure 9 Strain distribution on the capsule membrane around the capsule. The particles are numbered counterclockwise starting from the front end of the capsule.

In Figure 9, the distribution of strain around the capsule in all three channel geometries is presented. The particles are numbered from $i=1$ to $i=80$ in a counterclockwise manner starting from the front as indicated in Figure 8. The strain in all three cases is negative in the front and back and is positive at the top and bottom. The maximum positive strain in the capsule occurs in the $W=12r_c$ channel.

Figure 10 Shape evolution of capsule when deformed by microchannels of varying constriction widths (a) $W=8r_c$, (b) $W=12r_c$ and (c) $W=16r_c$.

In Figure 10, snapshots of the capsules are shown at different locations in the three channel geometries. In Figures 10(a), 10(b) and 10(c), the top row shows the capsule as it enters the constriction, the middle row shows the capsule in the middle of the constriction and the bottom row shows the capsule as it is exiting. Note that the snapshots are presented at different times since the capsules are transported rapidly in the widest channel and slowly in the narrowest channel. Due to effect of the channel geometry on the hydrodynamics, the capsule exhibits distinctly different shapes at the different locations. The capsule in the $W=8r_c$ channel forms a parachute shape from an ellipsoidal shape as it exits the constriction. The capsule squeezed through, $W=12r_c$ channel deforms a little less than $W=8r_c$ and forms a mushroom shape. The capsule squeezing through the $W=16r_c$ channel deforms the least amongst the three and exits with a bullet shape.

3.4 Pore Density in Capsule Membrane

It has been observed that the capsule membrane generally experiences a negative strain before entering the constriction. In the constriction, only parts of the capsule membrane near the channel wall experience a positive strain. A positive strain indicates the membrane is being stretched. As shown by TEM imaging of the membrane [38], disruptions in the membrane occur to accommodate the strain. These are temporary pores or SA channels. A recent study [33] on the formation of pores in Schlemm’s canal (SC) cells due to an applied equibiaxial strain suggests pore formation is a mechanosensitive process which is triggered by positive biomechanical strain. The number of pores formed is influenced by the magnitude of strain.

Figure 11 Time-averaged strain of the capsule in channels of three different constriction widths.

The average strain on the capsule wall is determined by averaging the positive strain around the capsule. The graph in Figure 11 shows the strain observed for different constriction widths as a function of capsule position in the channel. The capsule in the $W=8r_c$ channel has a lower average positive strain than the $W=12r_c$ channel. This implies that higher deformation may not necessarily result in larger regions of positive strain. Thus the $W=12r_c$ channel shows a higher averaged positive strain than $W=16r_c$ channel and the effect is maintained over a longer part of the channel in both cases. On the other hand, the effect of the constriction is felt on the capsule over a shorter length in $W=8r_c$ channel.

Figure 12 Piecewise linear fit of the pore density with strain to experimental data of [33].

As mentioned earlier, the formation of pores on the membrane is correlated to the positive strain experienced by the capsule membrane. In this work, we use this correlation to estimate the pore density in the capsule as it is transported in the channels. For this purpose we use the experimental data reported by [33]. We performed a piecewise linear fit to the data of the mean pore density of intracellular pores for 3 SC cell lines from human donors of different ages as a function of strain (Figure 12). For the strain range 0%–10%, the fit of the pore density as a function of strain is found to be ${\rho }_p(\epsilon )=5.75\epsilon +6.5$ and for strain in the range 10%–20%, the fit is ${\rho }_p(\epsilon )=0.4\epsilon +60$.

Figure 13 Time-averaged pore density on the capsule membrane in channels with different constriction widths.

Figure 14 Maximum strain on the capsule during deformation in three channel geometries.

We use this correlation to determine the pore density in the capsule as a function of position along the channel. The average ${\rho }_p$ for the three channel geometries $W=8r_c,12r_c,16r_c$ is given in Figure 13. The $W=12r_c$ channel has significantly higher ${\rho }_p$ due to positive strain in larger regions of the capsule membrane. The channels with $W=8r_c$ and $W=16r_c$ show lower ${\rho }_p$ as expected. This is in agreement with the experimental results [7], where the maximum mass transport was observed in constriction channel with $W$ as $50\%$ of the cell diameter.

Cell viability has been correlated to the maximum pore size developed in the membrane. Thus it should be ensured that with larger pore density, the maximum pore size should not be such that cell viability is affected. We plot the maximum strain experienced by the capsule as it is transported in the channels in Figure 14. If the maximum strain can be correlated to the maximum pore diameter, cell viability can be studied for designing better devices.

4 Conclusion

In this study, an FDPD model of the mechanoporation process in a microfluidic channel was developed. The strain in the capsule during transport in channels of different geometries was studied. It was found that the capsule experiences larger strain at intermediate channel widths. The strain persists over a longer length of the channels with wider constrictions. Next, a correlation of the strain with the pore density was used to predict pore density in the capsule during transport in the channel. This model can be used to study the effect of other parameters such as cell speed, time of deformation and other cell geometries to develop efficient devices for the diffusion of cargo into cells.

References

[1] M. P. Stewart, R. Langer, and K. F. Jensen, “Intracellular delivery by membrane disruption: Mechanisms, strategies, and concepts,” Chem. Rev., vol. 118, no. 16, pp. 7409–7531, 2018.

[2] A. Sharei et al., “Ex vivo cytosolic delivery of functional macromolecules to immune cells,” PLoS One, vol. 10, no. 4, pp. 1–12, 2015.

[3] G. L. Szeto et al., “Microfluidic squeezing for intracellular antigen loading in polyclonal B-cells as cellular vaccines,” Sci. Rep., vol. 5, no. February, pp. 1–13, 2015.

[4] J. Li et al., “Microfluidic-Enabled Intracellular Delivery of Membrane Impermeable Inhibitors to Study Target Engagement in Human Primary Cells,” ACS Chem. Biol., vol. 12, no. 12, pp. 2970–2974, 2017.

[5] A. Kollmannsperger et al., “Live-cell protein labelling with nanometre precision by cell squeezing,” Nat. Commun., vol. 7, p. 10372, 2016.

[6] T. DiTommaso et al., “Cell engineering with microfluidic squeezing preserves functionality of primary immune cells in vivo,” Proc. Natl. Acad. Sci., vol. 115, no. 46, pp. E10907–E10914, 2018.

[7] A. Sharei et al., “A vector-free microfluidic platform for intracellular delivery,” Proc. Natl. Acad. Sci., p. 201218705, 2013.

[8] J. Lee et al., “Nonendocytic delivery of functional engineered nanoparticles into the cytoplasm of live cells using a novel, high-throughput microfluidic device,” Nano Lett., vol. 12, no. 12, pp. 6322–6327, 2012.

[9] A. Sharei et al., “Plasma membrane recovery kinetics of a microfluidic intracellular delivery platform,” Integr. Biol., vol. 6, no. 4, pp. 470–475, 2014.

[10] M. Raab et al., “ESCRT III repairs nuclear envelope ruptures during cell migration to limit DNA damage and cell death,” Science (80-.)., vol. 352, no. 6283, pp. 359–362, 2016.

[11] T. Harada et al., “Nuclear lamin stiffness is a barrier to 3D migration, but softness can limit survival,” J. Cell Biol., vol. 204, no. 5, pp. 669–682, 2014.

[12] O. S. Pak, Y. N. Young, G. R. Marple, S. Veerapaneni, and H. A. Stone, “Gating of a mechanosensitive channel due to cellular flows,” Proc. Natl. Acad. Sci. U. S. A., vol. 112, no. 32, pp. 9822–9827, 2015.

[13] P. J. Hoogerbrugge and J. M. V. A. Koelman, “Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics,” EPL, vol. 19, no. 3. IOP Publishing, pp. 155–160, 01-Jun-1992.

[14] A. J. Ladd, “Numerical Simulations of Particulate Suspensions Via a Discretized Boltzmann Equation. Part 1. Theoretical Foundation,” J. Fluid Mech., vol. 271, pp. 285–309, 1994.

[15] D. H. Rothman and S. Zaleski, “Lattice-gas models of phase separation: Interfaces, phase transitions, and multiphase flow,” Rev. Mod. Phys., vol. 66, no. 4, pp. 1417–1479, 1994.

[16] W. Dzwinel, K. Boryczko, and D. A. Yuen, “A discrete-particle model of blood dynamics in capillary vessels,” J. Colloid Interface Sci., vol. 258, no. 1, pp. 163–173, Feb. 2003.

[17] M. M. Dupin, I. Halliday, C. M. Care, L. Alboul, and L. L. Munn, “Modeling the flow of dense suspensions of deformable particles in three dimensions,” Phys. Rev. E – Stat. Nonlinear, Soft Matter Phys., vol. 75, no. 6, p. 066707, Jun. 2007.

[18] H. Noguchi and G. Gompper, “Shape transitions of fluid vesicles and red blood cells in capillary flows,” Proc. Natl. Acad. Sci. U. S. A., vol. 102, no. 40, pp. 14159–14164, Oct. 2005.

[19] I. V. Pivkin and G. E. Karniadakis, “Accurate coarse-grained modeling of red blood cells,” Phys. Rev. Lett., vol. 101, no. 11, pp. 1–4, 2008.

[20] D. A. Fedosov, B. Caswell, and G. E. Karniadakis, “A multiscale red blood cell model with accurate mechanics, rheology, dynamics,” Biophys. J., vol. 98, no. 10, pp. 2215–2225, May 2010.

[21] W. Pan, B. Caswell, and G. E. Karniadakis, “A low-dimensional model for the red blood cell,” Soft Matter, vol. 6, no. 18, pp. 4366–4376, 2010.

[22] K. Müller, “In silico particle margination in blood flow,” Universität zu Köln, 2015.

[23] S. Z. Hoque, D. V. Anand, and B. S. V Patnaik, “The dynamics of a healthy and infected red blood cell in flow through constricted channels: A DPD simulation,” Int. J. Numer. Method. Biomed. Eng., vol. 34, no. 9, p. e3105, 2018.

[24] G. R. Lázaro, A. Hernández-Machado, and I. Pagonabarraga, “Rheology of red blood cells under flow in highly confined microchannels: I. Effect of elasticity,” Soft Matter, vol. 10, no. 37, pp. 7195–7206, Oct. 2014.

[25] X. Niu, T. W. Pan, and R. Glowinski, “The dynamics of inextensible capsules in shear flow under the effect of the natural state,” Biomech. Model. Mechanobiol., vol. 14, no. 4, pp. 865–876, Aug. 2015.

[26] B. Kaoui and J. Harting, “Two-dimensional lattice Boltzmann simulations of vesicles with viscosity contrast,” Rheol. Acta, vol. 55, no. 6, pp. 465–475, Jun. 2016.

[27] D. S. Hariprasad and T. W. Secomb, “Two-dimensional simulation of red blood cell motion near a wall under a lateral force,” Phys. Rev. E – Stat. Nonlinear, Soft Matter Phys., vol. 90, no. 5, p. 053014, Nov. 2014.

[28] X. Q. Hu, A. V. Salsac, and D. Barthès-Biesel, “Flow of a spherical capsule in a pore with circular or square cross-section,” J. Fluid Mech., vol. 705, pp. 176–194, 2012.

[29] S. Y. Park and P. Dimitrakopoulos, “Transient dynamics of an elastic capsule in a microfluidic constriction,” Soft Matter, vol. 9, no. 37, pp. 8844–8855, 2013.

[30] M. Y. Hwang, S. G. Kim, H. S. Lee, and S. J. Muller, “Elastic particle deformation in rectangular channel flow as a measure of particle stiffness,” Soft Matter, vol. 14, no. 2, pp. 216–227, 2018.

[31] K. K. Sreeja, J. H. Ipsen, and P. B. Sunil Kumar, “Monte Carlo simulations of fluid vesicles,” J. Phys. Condens. Matter, vol. 27, no. 27, 2015.

[32] B. Kaoui, G. Biros, and C. Misbah, “Why do red blood cells have asymmetric shapes even in a symmetric flow?,” Phys. Rev. Lett., vol. 103, no. 18, p. 188101, 2009.

[33] S. T. Braakman et al., “Biomechanical strain as a trigger for pore formation in Schlemm’s canal endothelial cells,” Exp. Eye Res., vol. 127, pp. 224–235, 2014.

[34] W. Pan, I. V. Pivkin, and G. E. Karniadakis, “Single-particle hydrodynamics in DPD: A new formulation,” EPL (Europhysics Lett., vol. 84, no. 1, p. 10012, 2008.

[35] A. Kumar, “Computational Modeling of Microfluidic Process using Dissipative Particle Dynamics,” University of Rhode Island, 2009.

[36] S. K. Ranjith, B. S. V. Patnaik, and S. Vedantam, “No-slip boundary condition in finite-size dissipative particle dynamics,” J. Comput. Phys., vol. 232, no. 1, pp. 174–188, 2013.

[37] A. Yazdani, M. Deng, B. Caswell, and G. E. Karniadakis, “Flow in complex domains simulated by Dissipative Particle Dynamics driven by geometry-specific body-forces,” J. Comput. Phys., vol. 305, pp. 906–920, Jan. 2016.

[38] A. R. Sharei, “Cell squeezing: a vector-free microfluidic platform for intracellular delivery of macromolecules,” Massachusetts Institute of Technology, 2013.

Biographies

Nishanthi Rajkamal received her bachelor’s degree in biotechnology in 2009, master’s degree in nanoscience and nanotechnology from Anna University in 2011. She is currently pursuing philosophy of doctorate degree in biomedical devices and technology at Indian institute of technology Madras. Her research is on modeling dynamics of mechanoporation in intracellular delivery devices.

Srikanth Vedantam is a professor in the Department of Engineering Design at the Indian Institute of Technology Madras. He received his ScD from the Massachusetts Institute of Technology in 2000. He has more than 20 years of teaching and research experience in the National University of Singapore and the Indian Institute of Technology Madras. His areas of expertise are in computational mechanics of solids and fluids, with a particular focus on discrete particle approaches.

European Journal of Computational Mechanics, Vol. 30_4–6, 409–430.
doi: 10.13052/ejcm2642-2085.30465
© 2021 River Publishers

Abstract

1 Introduction

2 DPD Model

2.1 Model of the Fluid

2.2 Model of the Capsule

2.3 Boundary Conditions

2.4 Model Parameters

3 Results and Discussion

3.1 Flow Through a Constriction Channel

3.2 Capsule Mechanics

3.3 Strain on Capsule Flowing in Channel with Varying Constriction

3.4 Pore Density in Capsule Membrane

4 Conclusion

References

Biographies