Modeling and Simulation of Multiphase Flow for Nanoparticle Translocation

: The delivery of bioactive compounds is often improved by their encapsulation within systems based on different materials, such as polymers and phospholipids. In this regard, one of the most attractive vehicles are liposomes, which can be produced by the self-assembly of phospholip-ids in aqueous buffered systems. Encapsulation of therapeutic magnetite nanoparticles (MNPs) within liposomes can be accomplished by direct translocation of their lipid bilayer by surface con-jugation of potent translocating peptides (and proteins) such as Buforin-II and OmpA. Here, we put forward the notion that to achieve reproducibility and optimize this process, it is possible to develop microfluidic systems that use flow-focusing methods to manipulate the interaction of suspended MNPs (ferrofluids) with the liposomes. With that in mind, we have developed an in silico approach to predict the performance of microfluidic devices specifically designed for the encapsulation process. This was done by running multiphysics simulations in COMSOL to evaluate the macroscopic flow of liposomes and suspended MNPs via a multiphase mixture model. Moreover, we estimated the corresponding interaction using a chemical reaction model based on embedding the Michaelis– Menten equation within the diluted species module’s transport. In this case, the enzymes-substrate interaction was considered similar to that of the MNPs-liposome. As a result, we were able to approach saturation kinetics that resemble that obtained experimentally for the uptake of functional-ized MNPs. Future work will be directed towards refining the model by considering more details on the possible stages during the interaction of the involved intermediates.


Introduction
Recent advances in nanoscale engineering have led to the development of cell-penetrating vehicles capable of internalizing different cell lines and even escaping endosomal routes of intracellular trafficking [1,2]. One of the most successful and promising nanoplatforms for therapeutics' delivery are magnetite nanoparticles (MNPs) mainly due to their biocompatibility, ease of synthesis and functionalization, low production cost, and marked responsiveness to magnetic fields [3][4][5]. To increase their properties as cell delivery agents, we immobilized the potent translocating protein OmpA as well as the peptide Buforin-II [6][7][8]. The obtained nanobioconjugates penetrated various cell lines quite effectively and showed exceptional endosomal escape levels. For this reason, we are considering their application in the delivery of therapeutic molecules for the treatment of neurodegenerative diseases and several cancer types. To achieve even higher cell penetration levels, avoid possible degradation processes, and extend life span and systemic circulation times, we decided to encapsulate the nanobioconjugates into liposomes.
Encapsulation of bioactive molecules into liposomes represents a potential approach for upholding their activity and stability under harsh environments [9]. Given their resemblance to cell membranes' phospholipidic composition, liposomes are considered one of the most attractive vehicles to study encapsulation processes and deliver therapeutics [10,11]. There are several methods to conduct encapsulation processes with liposomes, including mechanical dispersion, solvent dispersion, and detergent removal [12][13][14]. Most of these technologies are relatively mature; however, the obtained encapsulates exhibited large particle size distributions, varying morphologies, and, most importantly, led to relatively low encapsulation efficiencies [15,16]. A possible route to address these issues is to conduct the encapsulation with the aid of microfluidic devices. These technologies have proven useful to manipulate objects at the microscale by confining them to interactions with fluids within microchannels where the laminar flow regime is prevalent. Under such conditions, mixing processes can be controlled easily by altering the flow rates and the pitch and cross-sectional geometry of the microchannels. The emergence of microfluidics for controlled mixing and reaction processes has led to several applications in encapsulation, including drug delivery, high-throughput screening, and single cell analysis [17][18][19]. In the case of encapsulation into liposomes, recent attempts resulted in high throughput schemes; however, the efficiencies remain at relatively low levels [20].
With this opportunity for improvement, we started a research program dedicated to the design and manufacture of low-cost microfluidic devices for the encapsulation of MNPs into liposomes. We based our designs on the previous research work that demonstrates the ability of micromixers to manipulate suspensions of MNPs (i.e., ferrofluids) at a macroscopic scale. In this context, the ferrofluid mechanics and properties resemble those of an aqueous solution, allowing the use of these devices for particle focusing via fluid streamlines. Moreover, the role of channel geometry and external manipulation (such as using magnetic forces or acoustic fields) has been reported for altering those streamlines [21][22][23]. As for liposomes, their macroscopic properties can be modeled as those of water, following their assembly in an aqueous solution [24].
Here, we introduced an in silico approach based on a multiphase flow mixture model for estimating the level of interaction between the ferrofluid and liposomes. To model the encapsulation process, we implemented a model that couples the transport of diluted species and chemical reaction modules. This approach considers the Michaelis-Menten equation for the kinetics of interaction between liposomes and MNPs by referring to enzyme and substrate, respectively.

Mixture Model
To address the transport problem as a two-phase fluid, the liquid phase with the dispersed liposomes and the ferrofluid phase were modeled independently. Under this approach, each fluid phase has its own adjustable flow rate [25]. The relative flux of particles suspended in an aqueous medium is determined by the particle slip velocity [26]. Furthermore, the mixture model is a simplified model from the Eulerian-Eulerian model, which can consider n number of phases. This model assumes local equilibrium over short spatial length scales and that there are different velocity and concentration fields for each fluid phase involved. Besides, it also considers that there is a certain fraction of each phase within a given control volume [27]. This model solves the continuity and momentum equations for the mixture and the volume fraction equations for the secondary phases.
The continuity equation for the mixture model is given by (Equation (1)): where the subscript m refers to the mass-averaged variables of the mixture. The volume fraction of the particles and the liquid is defined as = ∑ / , (cell refers to a computational cell) and = 1 − , respectively.
Conversely, the momentum equation is obtained by summing up over the individual momentum equations for all phases, which is described by (Equation (2)): where , , and are the mass-averaged velocity, the mixture density, and the mixture viscosity, respectively (for all the involved phases). Finally, the transport equation corresponding to the solid-phase (particles) volume fraction is expressed as (Equation (3)):

Diluted Species Model
We approached the encapsulation problem by considering the interaction and subsequent encapsulation as the reaction between three different chemical species present in the mixture. Given that, this model also accommodates all types of materials transport that can occur through diffusion (modeled by the Fick's law) and convection, either alone or in combination with one another [28]. This model assumes that all species present are diluted in the solvent (i.e., the solvent's concentration is more than 90% mol) [28]. The model solves the conservation of chemical species, which is described by the following expression (Equation (4)): This equation considers both transport mechanisms (diffusion and convection). Where , ⃑ , and represent the concentration of the species i, the mass-average velocity vector, and the reaction rate expression for the species i, respectively. Lastly, is the diffusive flux vector given by (Equation (5)): where D denotes the diffusion coefficient and the concentration of the species i. In this work, we considered that the Ri term models the liposome-ferrofluid interaction through the Michaelis-Menten reaction rate expression. Thus, (Equation (6)) describes the reaction rate of each species present in the mixture: where Km is given by (Equation (7)): and Vmax is given by (Equation (8)): Considering this, R reflects the production of the liposome-particle complex following an irreversible reaction, resembling translocation of MNPs through the lipid bilayer of liposomes.

Computational Domain
The proposed computational domain corresponds to a simple microfluidic device consisting of two 2.12 mm diameter inlets and one 3 mm diameter outlet with a 30° intersection angle. Figure 1 shows the 2D computational domain.

Velocity and Pressure Profile
The device's velocity profile was calculated by implementing the Single-Phase Laminar Flow interface of COMSOL. Considering the device channel's dimensions, the inlets' flow rate = 30 mL/h, = 997 , = 8.9 × 10 −4 Pa·s for the water, we confirmed a laminar flow regime.

Mesh
Two meshes were considered for the entire geometry. The first, a Free Triangular mesh calibrated for fluid dynamics with a 0.118 mm and a 3.5 × 10 −4 mm for maximum and minimum element size, respectively. The second was calibrated for general physics with 8.3 × 10 −4 mm and 8.4 × 10 −7 mm for maximum and minimum element size, respectively. The small dimensions of both meshes were required to capture the dynamics of magnetite ( ≅ 10 nm).

Mixture Model Implementation
The Mixture Model, Laminar Flow (mm) interface from COMSOL was implemented to model the multiphase flow. The dispersed phase was considered as solid (Fe O ) magnetite nanoparticles with density = 5180 kg/m and diameter = 10 nm. No turbulence was considered, and a Schiller-Naumann Slip model was chosen. The continuous phase is considered as water with density = 997 kg/m and dynamic viscosity = 8.9 × 10 −4 Pa·s. For the mixture, a maximum packing concentration = 0.04 was specified. Accounting for the computational domain, a no-slip boundary condition was applied to the wall. The system's inlets were specified as the velocity flow rate = 30 mL/h while the dispersed phase inflow was set to = 0.02 volume fraction. For the outlet, a pressure = 0 boundary condition was specified to suppress backflow. In addition, an exterior dispersed phase condition of , = 0, as it was considered that all magnetite was contained within the system.

Diluted Species Model Implementation
Regarding the diluted species model, first, the laminar inflow and outflow conditions were specified as those used for Laminar Flow, with an inflow rate of 30 mL/h at each inlet, and an outflow driven by a pressure of 0. Then, liposome and particle inflow concentrations were set to the values listed in Table 1 through a parametric sweep. Further, liposome and nanoparticle diffusivity coefficients were assumed as 1.42 × 10 −9 m 2 /s and 7.19 × 10 −10 m 2 /s, respectively. Additionally, constant values and were set in the range of 0.05 to 1.2 and 0.05 to 3.0, respectively.
Moreover, considering the Michaelis-Menten kinetic equation (Equation (6)), a reaction condition was set by using the formation of an enzyme-substrate complex to model the liposome-MNPs interaction. Hence, each species' consumption rate was considered the negative of the product formation rate.

Results
Preliminary results of the velocity profile and the pressure field are shown in Figure 2.      As shown in Figure 4a,b, the liposome-MNP complexes are formed at about the middle of the main channel towards the liposomes' side. This profile largely superimposes with that of the MNPs (Figure 4c) but fails to do so with the liposomes (Figure 4d). In addition, as can be seen in Figure 4c, the MNPs' concentration profile remains largely unaffected, which can be explained by the fact that the concentration of liposomes limits the reaction. Figure 5 presents the liposomes uptake and the percentage of unencapsulated nanoparticles remaining in the fluid as a function of the initial MNPs' concentration and for different values of and . As shown in Figure 5a-c, by increasing the concentration of MNPs, the fraction of unencapsulated ones in the fluid increases steadily to about 60%, which is reflected in a decrease of the fraction of encapsulated particles. Moreover, Figure  5a-c suggest that by increasing or the saturation concentration, liposome uptake decreases. In contrast, the rate of uptake increases and Thus, and determine the uptake levels into liposomes that might be achievable for our nanobioconjugates.

Discussion
An initial approach for the interaction between MNPs (modeled here as a ferrofluid) and liposomes was achieved through a multiphase mixture model. The proposed model also allowed a more comprehensive understanding of the encapsulation process by considering a first-order reaction kinetics for the interaction of the liposomes and the MNPs. The obtained results suggest that the proposed microfluidic device fails to provide a complete mixture of both fluids (Figures 3 and 4c). This will negatively impact the extent of the encapsulation process and confirms the need for more intricate microchannel geometries where the interaction time is significantly extended. Because liposomes are the limiting reagent, the effectiveness of encapsulation is mainly driven by changes in their concentration profile behavior as corroborated by Figure 4a,b,d. Figure 5 shows that increasing the particles' affinity for the liposomes (through an increase in Km) fails to improve the uptake, as the fraction of unencapsulated MNPs remained at approximately 60%. This was also the case of the encapsulation rate, as an increase in Kcat was also ineffective. Taken together, these results also suggest that the model requires further refinement considering that it is possible to reach 100% uptake experimentally at higher MNP concentrations. A possible improvement is to include an intermediate step of interaction between the components.
Despite failing to reach the maximum uptake level observed experimentally, the model can predict how, by changing the affinity of the MNPs towards the liposomes, i.e., by changing Km, it is possible to obtain different saturation kinetics. The higher the affinity (or lower Km), the lower the concentration of MNPs required for saturation ( Figure 5), and the higher of MNPs remaining on the fluid. Similar results are obtained by increasing Kcat, which can be explained by faster translocation of the liposomes leading to higher uptake at low MNPs concentration.

Conclusions
Encapsulation of therapeutic MNPs into liposomes is a field that has gained considerable attention due to the possibility of improving even further their superior features as delivery vehicles. Here, we aimed at investigating whether this encapsulation process can be modeled effectively via multiphysics approaches. This was with the ultimate intention of designing and testing novel microfluidic devices in silico for high throughput encapsulation schemes. To accomplish this, we proposed a multiphase mixture model and a conservation of species model coupled to Michaelis-Menten kinetics. In the first case, we obtained mixing profiles that appeared insufficient for complete formation of liposome-MNP complexes, which was attributed to the simplicity of the contacting device that significantly limited the extent of interaction. This can be addressed in future work by devices with more complex geometries where also external sources of energy can be coupled to facilitate interactions (e.g., magnetic, acoustic, and thermal). The second model provided valuable information on the changes in the concentration profiles of the liposomes, MNPs, and the liposome-MNP complexes. In this regard, we confirmed the role of liposomes as limiting reagents in the proposed interactions. Additionally, it allowed us to investigate the role of affinity and reaction rate on the uptake of MNPs by the liposomes. Even though the model failed to reproduce the maximum level of uptake observed experimentally, it provided uptake curves that can be adjusted easily to experimental data by fine tuning of the Michaelis-Menten parameters Km and Kcat.