Experimental and Computational Study on the Microfluidic Control of Micellar Nanocarrier Properties

Microfluidic-based synthesis is a powerful technique to prepare well-defined homogenous nanoparticles (NPs). However, the mechanisms defining NP properties, especially size evolution in a microchannel, are not fully understood. Herein, microfluidic and bulk syntheses of riboflavin (RF)-targeted poly(lactic-co-glycolic acid)-poly(ethylene glycol) (PLGA-PEG-RF) micelles were evaluated experimentally and computationally. Using molecular dynamics (MD), a conventional “random” model for bulk self-assembly of PLGA-PEG-RF was simulated and a conceptual “interface” mechanism was proposed for the microfluidic self-assembly at an atomic scale. The simulation results were in agreement with the observed experimental outcomes. NPs produced by microfluidics were smaller than those prepared by the bulk method. The computational approach suggested that the size-determining factor in microfluidics is the boundary of solvents in the entrance region of the microchannel, explaining the size difference between the two experimental methods. Therefore, this computational approach can be a powerful tool to gain a deeper understanding and optimize NP synthesis.


INTRODUCTION
Polymeric nanoparticles (NPs) have attracted great attention in nanomedicine due to their favorable, tunable properties. 1−3 A wide range of drugs and diagnostic agents 4,5 can be conceptualized based on polymeric NPs. To increase their performance in drug delivery applications, researchers have attempted to control their size, shape, and surface charge to regulate the NPs' tissue penetration, detection, and clearance from the body. 6−8 A major aspect is the particle size, which arose as a crucial factor for in vivo applications and was determined to be ideal in the range of 10−200 nm. 9−11 However, control of the polymeric NP sizeespecially at the lower nanometer scaleis considered challenging. This can be overcome by sensitively regulating the interactions between polymers. Here, microchannels can be a good tool to tune the mass transfer between polymers. 12 In microfluidic platforms, the mixing step of polymer aggregation is done through molecular diffusion. 13 Many studies have been performed by researchers to tune the size, morphology, and structure of polymeric NPs via the design of microfluidic platforms. In this regard, various studies have shown the importance of microchannel design to control the size and size distribution of polymeric NPs. 14,15 Attempts have been made to optimize the synthesis parameters (e.g., flow rates for organic and aqueous phases) of poly(lactic-coglycolic acid)-poly(ethylene glycol) (PLGA-PEG) NPs in the microfluidic method. In vitro, those NPs with a defined structure and size have shown a higher uptake by cancer cells than macrophages. 16 Next to size, morphology of polymeric NPs has been engineered using the microfluidic approach. For example, Sadrabadi et al. 17,18 synthesized polymeric NPs using a microfluidic platform and highlighted its potential for achieving a spherical morphology. However, although the impact of the microfluidic synthesis approach on the size, morphology, and size distribution of polymeric NPs has been demonstrated, the underlying mechanisms are not fully understood.
Consequently, researchers have tried to study microchannels using computational fluid dynamics (CFDs). 19 The CFD studies have provided additional information about the distribution of concentration, velocity, and temperature in the microchannel that can be used to design microchips. 20 Experimental results and CFD calculations have provided good engineering information on the microfluidic effect on NP synthesis. However, how microfluidic effects take place at the molecular level is still unknown, and its atomic mechanisms pose many questions. What happens at the molecular level that leads to proper control of NP size in the microfluidic technique? What effect does the type of solvent molecules in the microchannel has on the particle size? What is the microscopic stability of NPs synthesized by the microfluidic method compared to other NPs? All of these questions remain unanswered.
Of course, investigations have been done to study microchannels using molecular simulation, but these works have not provided basic mechanisms, and most of the engineering and case studies have focused on the effect of microchannel applications in very specific cases. 21−23 Thus, there is still a gap of knowledge between the molecular mechanisms of the effect of microfluidic methods on the polymeric NP formation. Filling the gap could provide atomistic criteria for determining how microfluidics control the NP size.
Herein, we used a computational approach to understand size differences observed experimentally for bulk and microfluidic synthesis methods. In this regard, we synthesized small NPs via (i) a common bulk synthesis method and (ii) precipitated polymeric chains in microchannels (microfluidics). Analogously, in a computational approach, the micellization of NPs was simulated for bulk (random method) and microfluidic (interface method) syntheses. The simulation studies were conducted under the assumption that the synthesis of polymeric NPs in the microfluidic method is limited to the interface between solvents. The simulation results were matched with the experimental data to explain our observations atomistically. By presenting this computational approach, the current work helps better understand the effect of self-assembly of polymer strands and opens perspectives for improving microfluidic syntheses of polymeric NPs.

RESULTS AND DISCUSSION
To study empirical and computational aspects of the polymeric NPs' self-assembly in microfluidics, we selected PLGA-PEG polymers as amphiphilic systems that have been extensively investigated in microfluidic-based synthesis methods. 24 More-over, we considered the effect of riboflavin (RF) ligands on the polymeric NPs' dynamic and their impact on the final physicochemical properties. RF (Vitamin B2) and its derivatives have been widely used as targeting ligands in various studies in which they were shown to have high affinity, low toxicity, and efficient delivery both in vitro and in vivo. 25−27 The small RF-targeted NPs were synthesized by selfassembly of RF-conjugated PLGA-PEG chains with PLGA-PEG polymers in the bulk protocol, as reported previously. 28 Briefly, RF-conjugated polymer chains were synthesized first, and then, the proper ratio of the product was mixed with PLGA-PEG solution in acetonitrile (ACN) and added dropwise to water with gentle mixing ( Figure 1A).
In the current work, we aimed to model the self-assembly difference in bulk and microfluidic syntheses. For the latter, we decided to generate the NPs in a practical and easily accessible microfluidic device ( Figure 1B) with a square cross section of 20 μm in height (t-junction reactor). In this context, the same ratios of organic solutions were injected into microsized channels, while water was injected at the same time from both sides, all using syringe pumps. Similar methods have been widely utilized for the flow-based production of NPs by other groups. 12,29 Even though other more complex chip designs, such as herringbone mixers, allow for advanced engineering of particles compared to the t-junction reactor, 30,31 the MD method has only limited capacity related to the time and dimension of simulation boxes and is currently not able to adequately simulate more complex designs. Nevertheless, the formation of the NPs starts with the first contact between two phases, which can be simulated in a small (in our case, about 10 nm dimension) box. Within this small scale, only the interface of two phases would be considered, regardless of the chip design. Therefore, a rather simple experimental design was used that could still be mimicked in the computational approach.
Besides the hardware used for NP synthesis, different parameters can affect the NP size and configuration, such as concentration and flow rate. 32 For instance, Rhee et al. 33 Interface method used to model the self-assembly of RF-conjugated NPs in microfluidics. Blue and green molecules represent water and ACN, while RF-conjugated polymer chains are yellow, and PLGA-PEG strands are presented in gray.
reported an extensive analysis of NP size and size distribution based on the concentration of the polymer solution and the ratio of the organic solvent flow rate to the total flow rate (f). They concluded that with higher concentration and f values, NPs tend to aggregate and clog the microchannels. Moreover, at lower concentrations, NPs were shown to have smaller sizes. This is in line with our findings indicating that with lower concentrations (<5 mg/mL) NP sizes were very small. Unfortunately, clogging occurred, which might be explained by the fact that the NPs were small enough to penetrate the poly(dimethylsiloxane) (PDMS) structure. Therefore, swelling of the channel's wall and finally clogging happened. On the other hand, at higher concentrations (>15 mg/mL), aggregation and consequently clogging occurred ( Figure S1). The NP size increased with higher f ratios (>0.4). NPs with small sizes were reproducibly synthesized in 10 mg/mL concentration (organic phase) at f = 0.1, which were taken as the basic conditions for further studies ( Figure S2).
The computational approach to mimic the experimental conditions was based on a previous study 34 that used the MD method to simulate the self-assembly of targeted polymeric NPs considering the hydrodynamic behavior of solvents in microfluidic platforms. It has been reported that smaller sizes of NPs are produced in microfluidic setups due to the control over the reaction at the interfaces of parallel flowing streams. 35 In the microfluidic setup, we observed diffusive mixing of fluids only at the interfaces of parallel ACN and water streams in the microchannels, which gave rise to better control of precipitation of polymeric chains. In contrast, in the conventional bulk method, continuous mixing resulted in uncontrolled diffusion of water molecules to the ACN phase and precipitation of NPs. MD simulations were performed to explain the differences in nanoprecipitation in bulk and microfluidic syntheses.
(1) Random method: All molecules including water, organic solvent, and polymer chains were added randomly and uniformly into the simulation boxes ( Figure 1C). (2) Interface method: Water and ACN phases are formed where the polymer strands are located in the organic phase. Figure 1D represents the interface manner where due to the contact with water molecules at the interface, polymer chains aggregate into micelles. In both systems, the total number of polymer chains is 10, in which the number of PLGA-PEG chains' ratio to PLGA-PEG-RF chain numbers is defined as the PP/PPR ratio. Further mathematical inspection of the ligand impact was carried out by varying the PP/PPR ratios in both simulation approaches while keeping the total number of polymer chains constant. It should be noted that the PP/PPR ratio in the simulation is employed to represent the molecular weight content of the RFconjugated PLGA-PEG polymer in the NP composition. In total, 411 ACN molecules were used in both approaches, alongside 5000 water molecules. The copolymers used in the wet lab contain a huge number of monomers, which was not manageable in the molecular simulation. Consequently, simulated polymer chains had reduced monomer numbers compared to the actual synthesis but with the exact same molar ratio of lactide to glycolide as in the real condition.
where D is the diffusion coefficient of the solvent and R indicates the ratio of the flow rate of the polymeric stream to the total flow rate of water. Moreover, w f and w are the widths of the focused stream and channel, respectively. 12 According to the equation, the thickness of the focused stream is estimated to be ∼1.3 μm. As sketched in Figure 1D, two interfaces are formed at both sides of the focused stream. However, this exceeds the dimensional capability of the MD method. Therefore, we only considered one side of the interface in the simulation box. As the self-assembly of NPs is expected to occur analogously at both sides of the fluid interfaces, this simplification should not impair the comparison between experimental and simulation data. 2.1. NP Characterization in the Wet Lab. As expected, dynamic light scattering (DLS) results illustrate that all microfluidically synthesized NPs have a smaller size in comparison with those synthesized by the conventional bulk method (≤35 versus ≤45 nm). Figure 2A,B shows the NP size distribution for the nontargeted and 50 wt % RF-conjugated polymer samples, respectively. In both cases, NPs prepared by the microfluidic approach have smaller average hydrodynamic sizes and narrower size distributions. The transmission electron microscopy (TEM) image of nontargeted NPs from bulk synthesis ( Figure 2C(i, ii)) indicates that the particles have an average diameter of ∼26 nm (in the dry state), which is in good agreement with DLS results. Furthermore, the DLS results reveal that the NPs containing varying amounts of RFconjugated PLGA-PEG have an average size smaller than 50 nm. Interestingly, for both synthesis methods (microfluidic and bulk), the NP size drops to a minimum at 20 wt % of the RF conjugate. Further addition of RF-conjugated polymers in the composition results in larger sizes ( Figure 2D). Size differences between samples with a low content of RF in the compositions (bulk versus microfluidics) were statistically significant, while the size differences for 30−40 wt % RF were not significantly different. A more detailed inspection of each sample shows the largest size difference (∼13 nm) between both methods in the samples containing 10 wt % PLGA-PEG-RF, while the smallest (nonsignificant) difference was observed for the samples containing 30 wt % RF-conjugated polymers (∼3 nm). This might be explained by the higher molecular weight of the RFconjugated polymer chains and the assembling of the RF molecules on the surface. The latter might result in a different swelling behavior with a higher RF content that is less dependent on the synthesis method (see also Section 2.2.2).
To better understand the RF-dependent increase in hydrodynamic diameter by several nanometers, we examined the influence of RF on the surface of NPs. In accordance with other reports, 36,37 the ζ potentials of all NPs were negative (∼−27 to ∼−21 mV) ( Figure 2E). Surface charges of the nontargeted NPs from either process are approx. −23 and −24 mV for bulk and microfluidic syntheses, respectively. Interestingly, for conventionally synthesized NPs in the presence of RF, the ζ potential shows a slightly different value (only a 3 mV increase) and the tendency toward neutral charges. This could also affect the stability of the NPs and consequently their pharmacokinetics and pharmacodynamics in a biological environment. Even though not statistically significant, the difference in the ζ potential encouraged us to conduct further mathematical analysis in the atomic scale to explore the invisible features of the structures.
2.2. Molecular Dynamics. 2.2.1. Radius of Gyration. The radius of gyration (R g ) is the average distance of molecules from their center of mass (eq S2). 28,38,39 R g shows the radius of the NPs formed in the simulation box and has a concept similar to the hydrodynamic radius (R H ). 40 With increasing attraction energy between the molecules, they move closer to each other. Thus, Rg can represent the attraction energies between the molecules.
To explain the results, we present a more detailed information on the R g . As mentioned before, during the assembly of the polymer structure, polymer strands change their locations. Therefore, R g of the NPs changes during the simulation to reach a stable condition in the simulation and to form an NP. Figure 3A displays the variation of R g for PP/PPR ratios of 8:2 for the random (bulk) and interface systems. It reveals a steady decrease in the first 10 ns, while the diagram increases sharply until t = 50 ns and then keeps slowly declining with passing time. Strikingly, the trend observed for both random and interface systems is similar, with a small difference in the R g at each nanosecond. However, in the last 5 ns, fluctuations of the interface system are smaller than in the random simulation. There is an obvious difference between the two methodologies that had the same interacting number of polymer chains. We attribute this difference to the presence of water and ACN molecules between chains during the precipitation, and consequently a longer distance between chains, which resulted in a larger R g for the random approach.
Similar diagrams are plotted and presented in Figure S3 for other PP/PPR ratios.
To provide comparable results with the empirical section of the paper, we present the R g values for samples in the form of the normal distribution of R g reported during the progress of the simulation. Figure 3B,C is chosen to represent the shift in the R g when simulating the same polymer composition in the random design or interface mode. Figure 3D represents the final R g values that were achieved for the NPs simulated within random or interface precipitations. The results also reveal favorably smaller R g for the interface approach, which is shown with the solid line, i.e., molecular diffusion in the interface of the organic phase with water plays an important role in controlling the NP diameter and properties through the flow-based synthesis. In a previous report, 28 we used a simplified system in which the polymer strands (in the absence of organic solution) come in contact with the water molecules and formed NPs. In this article, in line with other reports, 41,42 the effect of the organic solvent was neglected in the simulation, and the simulation was performed for the interaction of polymers with water molecules. Interestingly, in our simulation, the random manner leads to similar results. Herein, to have comparable preconditions to the wet lab, in the interface as well as random approaches, ACN molecules were added into the simulation box. Interestingly, computational approaches with and without ACN generated similar R g values.
Overall, stable and compact compositions occurred at the PP/PPR ratio of 8:2 for both bulk and interface modes. We hypothesize that the self-assembly of polymeric strands at the interface of the solvent and antisolvent controls the structure of the products, which is regarded as a phenomenon that regulates the mass and heat transfer in microfluidics. The visual examination of the particles with an 8:2 PP/PPR ratio confirms the more compact and tidy NPs resulting from the interface approach ( Figure 3E). Consistent with Wilkosz et al., 41 particles contain no water molecules inside the mass, neither in the interface mode nor in the random method, which is attributed to the formation of hydrogen bonds between ester groups and water molecules.
2.2.2. Properties of Simulated Nanoparticles. To further evaluate the outcome of the computational approach, the simulated stable NPs were visualized without surrounding water molecules (Figure 4). Both particles contain ACN molecules in their mass, which were entrapped within the polymer chains due to the interaction of PLGA polymers with the organic solvent. It should be mentioned that, as we reported before, 34 RF-molecules were located at the surface of particles and both assemblies (interface and random methods) contain no water in their mass. The amphiphilic polymer coil in the NP forms a curved shape; the PLGA head is in the core and the RF tail bends toward the center. This can contribute to the hydrophobic interactions of RF with PLGA. Nonetheless, RF is located at the surface of the NPs, which is attributed to the steric hindrance because of the neighbored PEG chains.
To understand the compactness and tidiness of particles generated at the interface, we examined the density of polymers and ACN along with the simulation box ( Figure  4B) with respect to Figure 4A. As can be seen, the density of both increased inside the simulation box at the interface of water/ACN, where the precipitation takes place. However, in the case of bulk design, the distribution of the polymer chains is bimodal with the maximum amount at the 2 and 4.5 nm coordination, while the peak of ACN density occurs at approx. 3.1 nm. These observations confirm the extensive config- uration of the particle formed in a random manner with larger R g values.
To provide a quantitative description of the study regarding the impact of the interface approach on the compactness of NPs, we explored the solvent-accessible surface area (SASA) of the NPs in the simulation. For instance, at an 8:2 ratio, the SASA measurements for simulations of the interface and random NP formation provided values of 272 and 256 nm 2 at 100 ns, respectively. Like other parameters, SASA changes in the course of the simulation. We studied the SASA for all samples; Figure 4C presents two cases (NPs containing two and four RF-conjugated polymers) as examples. As indicated by this figure, during the precipitation, the difference between the two cases (8:2 and 6:4 ratios) in the random mode is bigger than the difference in the same samples in the interface style. It is obvious that for the interface samples, hydrophobic segments of polymer strands are more buried in the mass, and therefore, less surface area is available for solvents. In contrast, for the random samples, polymer strands are more scattered and exposed toward the solvent, resulting in higher SASA values.

Energy Analysis.
To further explore the simulation results, energy analysis was performed on all systems. The Gibbs free energy is the absolute criterion of stability. In the self-assembly of NPs (including polymeric NPs), the vdW attractions and electrostatic repulsions are the main components of the effective energy interactions. 43,44 For the presented system, the main parts of the Gibbs free energy are electrostatic and vdW energies. As higher absorption interactions induce stronger aggregation, numerous scientists have studied these energies to predict the attraction and repulsion between molecules and their stability. 45,46 Figure 5A shows the results of energy analysis in simulations considering random and interface syntheses with a PP/PPR ratio of 8:2. Due to the neutral charges of the polymer strands, for both approaches (random and interface), electrostatic interactions have smaller contributions to the total energy. Therefore, the vdW attraction is the dominant energy type in the NP formation. For instance, Figure 5A(i) shows that in the random approach, the electrostatic interactions start at −190 (kJ/mol) and end at −980 (kJ/mol), while vdW attractions continuously drop from −970 (kJ/mol) to −3300 (kJ/mol). The slopes of vdW and total energy are steep in the first 20 ns. After that point, although there are fluctuations, an equilibrium state is reached at the end of the simulation. The diagram confirms that NPs assemble within the 20 ns, and the fluctuations are due to the conformational and locational changes of the polymer strands in the NPs at the later phase of the simulation.
We explored cases of energy contribution at the end of the simulation ( Figure 5B). It is noteworthy that in most cases, the vdW attractionand consequently the total energyis stronger for interface systems compared with those of the random ones, i.e., being the main reason leading to the smaller R g values in the interface approach. Consequently, this explains why NPs produced in microfluidics have more compact structures. 47 Furthermore, energy interaction reaches the minimum at an 8:2 ratio (the most stable samples) in either the interface or random style of simulation.
We carried out a similar analysis for other cases and noticed that for all cases (either interface or random for all PP/PPR ratios), the contribution of electrostatic repulsions is about 20% of the total energy calculations. This means the vdW attraction is the driving force of the NP formation, preventing the electrostatic repulsions of the disassembled NPs. With further increase in the presence of RF-conjugated polymers, the vdW contribution is decreased and electrostatic repulsions become dominant, which indicates that the RF-conjugated polymers repel each other (see the Supplementary Information Figure S4).
Our previous study 24 led us to wonder about the role of solvents in the self-assembly of NPs in microchannels, which can also serve as another piece of evidence on the reliability of the pathway. Consequently, we repeated one of the interface simulations replacing ACN with dichloromethane (DCM). DCM is hydrophobic, switching the self-assembly process to an oil/water emulsion for the microfluidic approach. 48 As expected, this led to the formation of a layer of polymer chains that were attached to the interface of DCM and water without forming any NPs in the simulation ( Figure 5C). As stated previously, the solvent has a decisive role in microfluidic synthesis. In other words, replacing a water-soluble solvent with an insoluble solvent can shift particles' size from nano-to micro-scale, 24 since it changes the formation process. Energy analysis of the DCM/water interface underlined this assumption as interaction energies among polymer chains drastically increased. For instance, the total energy at the end of the simulation course is about −100 kJ/mole ( Figure 5D), while in a similar case with the ACN solvent, it is approx. −3000 kJ/mole. Therefore, polymers can neither attract each other nor form NPs and are collected as a layer of polymers in the interface.
Gibbs' free energy calculations (as a criterion for determining the stability of these aggregations) were performed using the umbrella sampling technique. Figure 5E demonstrates the results for both simulation approaches, where samples at the PP/PPR ratio 8:2 had the minimum Gibbs free energy. It can be concluded that the higher the molecular attraction energy between the polymers, the more stable the NPs. The results are in line with previous experimental reports, 47 indicating that synthesis in the microchannels improves particles' stability, fosters stable polymeric aggregation, and produces NPs with a more compact morphology and smaller sizes. 49,50 We could confirm these results by showing that the particles produced in the microfluidic system had stronger intermolecular attraction forces resulting in smaller polymeric NPs. The strong intermolecular attraction both stabilizes and determines the size of the NPs.

Mean-Squared Displacement.
In the search for a logical explanation for the observed results (experimentally and computationally) and differences, we used mean-squared displacement (MSD) to study the molecular mobility of polymeric chains during the self-assembly. The slope of the MSD during the process is proportional to the diffusion coefficient of molecules. 51,52 Therefore, a larger slope of the MSD diagram indicates a longer displacement and higher molecular mobility of the considered molecules. As shown in Figure 5F, employing the random approach in the assembly of similar compositions (8:2 ratio of PP/PPR strands) leads to a substantial displacement of polymer strands in comparison with the interface approach. On the other hand, the interface of the organic solution/water strongly controls chains' immigration during the assembly. The MSD scrutiny of NPs at all PP/PPR ratios ( Figure S5A) suggests a better control for the molecular diffusion of polymeric strands in the interface method, resulting in shorter distances between the assembled polymer chains. Consequently, the formed NPs were tighter, and molecules attracted each other strongly. Moreover, the MSD analysis revealed that the method of simulation (random and interface) had no significant impact on the molecular diffusion of water or the organic solvent ( Figure S5B,C), suggesting that the diffusion of polymeric chains at the interface controls the particle formation.
2.3. Comparison between Experimental and Computational Results. The trends observed by R g analyses in both simulation approaches are in agreement with the mean NP sizes determined in experiments using DLS measurements ( Figure 6). The minimum sizes of the hydrodynamic radii (R H ) of synthesis approaches (microfluidic and bulk) occurred at the same PP/PPR ratios as in the computational approach. Nevertheless, there are a few differences between experimental and computational analyses: for instance, in the computational analysis, nontargeted particles had very similar R g values when simulated in the interface mode and the random method, while in experimental synthesis, the same composition of polymers led to NPs with approximately 10 nm difference (hydrodynamic radii as measured by DLS). We attribute the difference to our simplified conditions in the simulations, as well as to the unknown interactions among the chains or aggregation during the synthesis procedure.

CONCLUSIONS
In conclusion, we prepared RF-conjugated polymeric NPs by a conventional bulk method and a microfluidic platform. Samples with the same molecular composition yielded smaller sizes in microfluidic synthesis compared to bulk synthesis. Both synthesis methods resulted in particles with the smallest hydrodynamic diameter at 20 wt % of the RF-conjugated PLGA-PEG in the composition of the NPs. The differences between microfluidics and conventional self-assembly of nanoscaled polymeric particles were further evaluated using classical molecular dynamics simulation. In this regard, two different modes of simulations (random and interface) were designed corresponding to the bulk and microfluidic syntheses in the wet lab, respectively. Initial evaluations revealed that at the same molecular composition, samples assembled by the interface method have a smaller R g in comparison with aggregates generated by the random method. Moreover, the computational results indicated that for both methods, the most compact and stable particles were generated at an 8:2 PP/PPR ratio corresponding to 20 wt % NPs with RFconjugated PLGA-PEG. It has been observed that vdW interactions mostly contributed to the assembly process. Scrutinizing mean-squared displacement of polymeric chains revealed that molecular mobility at the interface method is controlled, which is the main reason for the controlled diffusion of the polymers that give rise to more compact particles when compared with the random mode.
The proposed interface mechanism can be used for atomistic simulation, analysis, and prediction of other flow-based synthesis methods. Moreover, the presented results provide a logical basis to simplify simulations. Thus, this manuscript contributes to our molecular understanding of the selfassembly phenomenon, which is invisible and inaccessible via experiments for conventional and microfluidic methods. 4.2. Synthesis of Riboflavin Derivative. Riboflavinconjugated polymers were synthesized according to a modified procedure reported by Tsvetkova et al. 25 Briefly, derivation of a proper functional group on the riboflavin molecule was carried out as described next (Supporting Information, Scheme S1).
First, 6.0 g of riboflavin (∼16 mmol) was dissolved in 80 mL of acetic acid and cooled to 0°C. Perchloric acid (200 μl) was added to the reaction mixture, followed by the addition of 40 ml of acetic anhydride, and the temperature was increased up to 100°C and stirred overnight. The mixture was diluted with water and extracted with chloroform. The organic phase was dried over anhydrous sodium sulfate and filtered. The solvent was removed under a vacuum. Then, 8.5 g of 2′, 3′, 4′, 5′-tetra-O-acetylriboflavin was obtained as a yellow product (98% yield). This 8.5 g of 2′, 3′, 4′, 5′-tetra-O-acetylriboflavin (15.6 mmol, 1 equiv) was dissolved in anhydrous acetone. Then, 4.35 ml of ethyl bromoacetate (39 mmol, 2.5 equiv) and 8.63 g of potassium carbonate (62.4 mmol, 4 equiv) were added to the mixture and stirred overnight at room temperature. Subsequently, the product was diluted with 0.1 M HCl and extracted with dichloromethane.
The organic phase was removed, as described in the previous step. The achieved substance was purified by column chromatography on silica gel using DCM/MeOH (gradient from 0% MeOH to 10% MeOH) as the eluent. Subsequently, the organic solvent was removed under a vacuum.
The obtained product was heated to reflux in 2 M HCl for 2 h at 100°C. The solvent was removed under reduced pressure to obtain carboxymethylriboflavin.
4.3. Synthesis of PLGA-PEG-RF. The riboflavin derivative was conjugated to the amine group of PLGA-PEG-NH 2 polymer chains using EDC coupling. PLGA-PEG-NH 2 (0.01 mmol, 1 equiv) and carboxymethylriboflavin (52 mg, 10 equiv) were dissolved in DMSO, followed by the addition of EDC (23 Figure 6. Comparison of the empirical and computational results, which highlights the accuracy and reliability of both simulation approaches. mg, 10 equiv). The reaction mixture was stirred overnight at room temperature. Finally, after dialysis using the membrane tubing (MWCO 3 kDa) and from freeze-drying, a crude yellow polymer powder was obtained. To ensure purification, the crude product was dissolved in the acetonitrile and filtered using syringe filters with a 0.20 micron PTFE membrane. The reaction yields were measured as 76 mg, 76%.
4.4. Microfluidic Chip Fabrication. For the preparation of a microfluidic chip, the Sylgard 184 silicone elastomer poly(dimethylsiloxane) (PDMS) was mixed with the included curing agent the ratio of 10:1 and mixed vigorously for 5 min. The achieved mixture was cast on the master and placed in a desiccator to eliminate entrapped air. To finalize the replica molding process, the PDMS was cured in an oven for 2 h at 60°C . The obtained replica was cut and peeled from the master, and holes were punched to provide a connection to the outer world. The PDMS slab and glass slide were cleaned via sonication (10 min), dried, and subsequently exposed to an oxygen plasma (30 s, 37 W, 40 mL/min) and united. Post bake of 1 h at 60°C was performed to increase the bonding strength. 53 4.5. NP Preparation. PLGA-PEG-RF NPs were produced by the nanoprecipitation method in conventional bulk and microfluidic environments. PLGA-PEG and PLGA-PEG-RF were dissolved in acetonitrile at different ratios with a final concentration of 10 mg/mL. 4.5.1. Bulk Synthesis. For each sample, 400 μL of the polymer solution was diluted in Milli-Q water dropwise to 1 mg/mL, and the mixture was stirred for 2 h. The residue of the organic solvent was evaporated overnight at room temperature.
4.5.2. Microfluidic Synthesis. For each formulation, 400 μL of the polymer solution in acetonitrile was introduced into the chip as core flow and focused from both sides with Milli-Q water. The flow rates for the polymer solution and water were adjusted to 500 and 4500 μL/h, respectively. In all experiments, polymer solutions were diluted to the final concentration of 1 mg/mL. The residue of the organic solvent was evaporated overnight at room temperature.
In a first pilot study, we synthesized NPs with various concentrations and with different flow rate ratios of the organic phase to the total flow rate. These initial analyses were carried out using nontargeted polymers (PLGA-PEG). At lower and higher concentrations, swelling of the channel's wall and consequently clogging of the microchannels happened ( Figure  S1). Figure S2 represents the size and polydispersity index (PDI) of the prepared NPs. With an increase in the f ratio, size as well as PDI increases. At a polymeric concentration of 1 mg/mL, the synthesis procedure was stopped in most f ratios due to clogging. Therefore, the related data are not reported in the diagram. With an increase in the concentration, at the same f ratio, NP sizes increase.
4.6. Characterization. 4.6.1. Transmission Electron Microscopy (TEM). The TEM images of the polymeric NPs were acquired by employing a CM200-FEG-Philips TEM device (FEI Co., Eindhoven, The Netherlands). To this end, a dilute suspension of NPs was prepared and deposited onto the Cu grid with a carbon film. The morphology and size of the particles were characterized via diffraction (amplitude) contrast and (for crystalline materials) through high-resolution (phase contrast) imaging. 54 TEM images were analyzed using ImageJ 1.52v software.
4.6.2. Dynamic Light Scattering (DLS). Dynamic light scattering was performed at a fixed angle of 173°using a Zetasizer (Zetasizer 3000HS, Malvern Instruments Ltd., Worcestershire, UK) for the diluted suspensions in Milli-Q water at 25°C. For each sample, 100 μL of the sample was diluted with water in a disposable cuvette, and the size and ζ potential of NPs were measured.
GraphPad Prism 9.1.2 was used to perform an unpaired Student t-test. The results were considered significant if P < 0.05. 4.7. Molecular Dynamics Simulations. In the simulation part of this work, at the inlet region of the microchannel, due to the shear stress of the flow and the surface tension of the solvents, an interface is formed between acetonitrile and water. Polymers are initially soluble in acetonitrile and diffuse to the interface between acetonitrile and water due to the concentration gradient. Thus, the aggregation of polymers occurs at the interface, and the interface controls the mass transfer and formation of polymeric NPs. As a result of the complete mixing of solvents and the increase in molecular collisions, it is assumed that the size of the NPs does not change after the microchannel inlet area. Therefore, in general, we have proposed a mechanism for this system, called the interface, according to which the interface simulation of solvents can be a good criterion for simulating the whole microchannel.
Molecular dynamic simulation was performed to simulate the nanoprecipitation process. The used copolymers contain a large number of monomers, which was not manageable in the molecular simulation. Consequently, simulated polymer chains had less monomers than real (actually used) polymers: the ratio of PLGA/PEG monomers in the PLGA 5kDa -PEG 3.4kDa -Rf polymer was 38:77 in the experimental mode and 17:34 in the simulation mode. The ratio for the PLGA 3kDa -PEG 2kDa polymer is, however, almost identical, with 23:45 in the experimental mode and 10:20 in the simulation mode. Additionally, PLGA polymers had an equal molar ratio of lactide to glycolide in the simulation as in the real condition. Figure S6 shows the structure of monomers and constructed riboflavin, polymer molecules. Table S1 shows the charge and type of atoms in each molecule. The polymer chains had 26 and 30 nm length for PLGA 3kDa -PEG 2kDa and PLGA 5kDa -PEG 3.4kDa -Rf, respectively. Here, the chains were sequentially optimized using Avogadro and HyperChem8.0 software. 55 Next, optimization steps were performed using CP2K software (b3lyp and basis set of 3-21G+). At the main optimization step, we used GROMACS 2019.5 (OPLS-AA force field31) in the EM (10 kJ/mol/nm minimum force), NVT (constant number of atoms, N; constant volume, V; constant temperature, T) (500 ps), NPT (constant number of atoms, N; constant pressure, P; constant temperature, T) (500 ps) and then MD (100 ns) simulation in 2 fs time steps. All cases contained 30 000 water molecules (SPC/E water model) in boxes with a dimension of 3 × 3 × 30 nm 3 .
The basic force field includes bonded and nonbonded energies (E total = E bonded + E nonbonded ) where each has its components 56 Therefore, it can be stated that the total enthalpy changes for this system (in which there is no bonded interaction in the selfassembly of polymers) are equivalent to the sum of the changes of the van der Waals (vdW), electrostatic, and hydrogen bonds. The vdW energy is calculated based on the Lennard−Jones equation. Here, E is the intermolecular potential between the two atoms or molecules. σ is the distance at which the intermolecular potential between the two particles is zero. σ gives a measurement of how close two nonbonding particles can get and is thus referred to as the van der Waals radius. It is equal to one-half of the internuclear distance between nonbonding particles. r is the distance of separation between both particles (measured from the center of one particle to the center of the other particle).
The energy of electrostatic interactions follows Coulomb's law F is the electric force between the two charges q 1 and q 2 and r is the distance of separation. K is the Coulomb constant and ε 0 is the permittivity of space. The cutoff radius was adjusted at 1.4 nm for the (vdW) and Coulomb interactions. We used the Coulomb energy algorithm and the particle mesh Ewald (PME) method. The MD simulation was carried out with the isotropic Parrinello− Rahman algorithm 56,57 at 1 bar and with the Nose−Hoover method (velocity-scaling algorithm in NVT and NPT) at 300 K. The constraint algorithm was based on the Lincs algorithm, which is only used for the hydrogen bonds. The partial atomic charges of the structures are calculated using CP2K software 58 (pop = esp). Finally, for the main simulation, we do the same option of EM, NVT, NPT, and MD as before.
4.8. Gibbs' Free Energy Calculation. The Gibbs free energy (ΔG = ΔH−TΔS) was calculated for all samples by the umbrella sampling technique, as described previously. 28,59 Briefly, aggregated MD simulation structures were used for the input of umbrella sampling simulations. The computational analysis included two simulation steps. First, the pull code was used to separate one of the polymers from the NP. Subsequently, 100 configurations (each as 1 nanosecond) were extracted from the pull code simulation. After applying the pull code for the polymer strand, it restrained at increasing the center-of-mass (COM) distance from polymer strands, which led to the generation of various configurations for each location. The probability mass function (PMF) curve was extracted in the restrain stage using the polymers strands' positions to the COM. In other words, integration of PMF was performed corresponding to the series of configurations. Finally, the Gibbs free energy was obtained by the weighted histogram analysis method (WHAM) 60,61 on all configurations. The WHAM analysis method is a very powerful technique based on the estimation of the statistical uncertainty of the probability distribution provided by the umbrella method. 62,63 As a result, the smallest uncertainty can be computed by PMF results. The dimensions of the simulation box were 10 × 10 × 30. The energy minimization was applied to the NPs, and then, the NPT equilibrium was performed on the molecules. The temperature and pressure algorithms were the Nose−Hoover and Parrinello−Rahman, respectively. In general, in umbrella simulation, a polymer was separated from the aggregated polymer set by the pull code method. In total, 100 configurations were extracted from the separation steps of this polymer. Then, an MD simulation was performed on every 100 configurations. Finally, the WHAM algorithm was calculated on all configurations and the Gibbs free energy was obtained.
Micrographs of microchannels in the pilot study ( Figure  S1); synthesis of NPs with various flow rate ratios ( Figure S2); radius of gyration for structures with variable amounts of RF-conjugated polymers ( Figure  S3); energy contribution in all simulated cases ( Figure  S4); mean square displacement (MSD) analyses ( Figure  S5); synthesis and characterization: three-step synthesis route of carboxymethylriboflavin including NMR spectra for all synthesis products; computational approach:molecular dynamics simulation; synthesis and characterization: gyration radius: RMSD analysis; illustration of molecules used in simulations ( Figure S6); specification of the OPLS-AA atom types and partial charges used in simulations (Table S1); comparison of the obtained results with the repeated simulation (