Fluid Dynamics Optimization of Microfluidic Diffusion Systems for Assessment of Transdermal Drug Delivery: An Experimental and Simulation Study

: Organ-on-a-chip technologies show exponential growth driven by the need to reduce the number of experimental animals and develop physiologically relevant human models for testing drugs. In vitro, microfluidic devices should be carefully designed and fabricated to provide reliable tools for modeling physiological or pathological conditions and assessing, for example, drug delivery through biological barriers. The aim of the current study was to optimize the utilization of three existing skin-on-a-chip microfluidic diffusion chambers with various designs. For this, different perfusion flow rates were compared using cellulose acetate membrane, polyester membrane, excised rat skin, and acellular alginate scaffold in the chips. These diffusion platforms were integrated into a single-channel microfluidic diffusion chamber, a multi-channel chamber, and the LiveBox2 system. The experimental results revealed that the 40 µ L/min flow rate resulted in the highest diffusion of the hydrophilic model formulation (2% caffeine cream) in each system. The single-channel setup was used for further analysis by computational fluid dynamics simulation. The visualization of shear stress and fluid velocity within the microchannel and the presentation of caffeine progression with the perfusion fluid were consistent with the measured data. These findings contribute to the development and effective application of microfluidic systems for penetration testing.


Introduction
For testing drug diffusion through physiological barriers, different instruments are available.For topical/transdermal drug delivery, the vertical and horizontal static Franz diffusion cells are widely used.In the last decade, however, many dynamic systems have also been developed.In these microfluidic devices, microchannels are fabricated and continuous perfusion is provided to mimic the fluid flow in the living tissue.This fluid flow results from capillary microcirculation, lymph microcirculation, and the movement of extracellular fluids.In the traditional method of microdialysis, diffusion also occurs in the target tissue, where the microdialysis probes are implanted [1].A molecular exchange takes place through the semipermeable membrane of the probe in the extracellular matrix in the direction of the concentration gradient.The recovery of the system for a certain molecule depends on several factors, including the flow rate of the perfusion fluid, which can be an artificial cerebrospinal fluid or artificial peripheral perfusion fluid (both fluids are water-based salt solutions); the temperature; and the diffusion surface.A lower flow rate results in higher recovery of the test compounds and increasing the flow rate leads to an exponential reduction in the probe relative recovery.Therefore, in an optimal microdialysis setup, a 0.5-2.0µL/min flow rate [2,3] is applied to reach the maximum solute concentrations in the dialysate samples and leave enough time for the molecular dynamic movements and the diffusion of the test compound.On the other hand, under physiological conditions, the blood flow rate is the highest in the arteries (3-26 mL/min, depending on the diameter of the blood vessels), followed by the veins (1.2-4.8 mL/min), and is the slowest in microcirculation, where the diameter of the blood microvessels in humans is 5-10 µm [4].The capillaries present the highest cross-sectional surface for fluid perfusion in the body.
The currently available organ-on-a-chip or skin-on-a-chip systems apply different perfusion flow rates depending on the microchannel design and the placement and integration of cells or ex vivo tissues in the devices [5][6][7][8][9].To achieve physiologically relevant conditions, the shear stress on the channel wall and cells and the tissue preparations should be considered.In organ-on-a-chip microfluidic systems, cell cultures or human tissues like skin can be integrated into the device under physiological conditions.The systems can be designed with different complexities using one or more cell types and mimicking physiological or pathological conditions (e.g., accelerated aging, etc.).Also, the size of these micronized testing platforms makes it possible to reduce the amount of consumables, cells, tissues, test articles, etc.The relatively low cost and ease of manufacturing make this technology attractive for research purposes.The key considerations in designing in vitro microfluidic devices for drug testing are the (1) selection of cell types or tissue preparations for testing, (2) optimization of the design of the microfluidic channels, (3) optimization of the raw materials selected for the construction of the microfluid device, (4) selection of the best manufacturing technology, and (5) adjustment of the conditions like temperature, flow rate, shear forces on the cells or tissues, medium composition, oxygen content, etc., to produce physiological conditions.
In the current study, three different devices were compared and tested from a fluid dynamics point of view.Two skin-on-a-chip devices (single-channel and triple-channel setups) were developed in-house by our institution, and one was purchased from IVTech Srl (LiveBox2) (Ospedaletto, Italy) [10][11][12].To design better testing platforms and optimize those currently used, various flow rates were experimentally investigated, and mathematical simulations and computational fluid dynamic calculations were performed.The aim of this study was to determine the best testing conditions for topical drug delivery in various diffusion chambers.To evaluate the degree of diffusion of a 2% caffeine-containing cream, a hydrophilic model formulation was utilized on the surface of the synthetic membranes (cellulose-acetate and polyester), excised rat skins, and 3% alginate scaffolds.The alginate will be used in later experiments as a bioink scaffold for the 3D bioprinting of artificial human epidermis.The current series of experiments were conducted to optimize the conditions for future human skin 3D bioprinting studies.

Caffeine Cream Formulation
The caffeine cream was prepared with the following composition (Table 1): The cream was prepared ex tempore under magistral conditions.The lipophilic components were melted in an enamel bowl over a water bath, polysorbate 60 was added, homogenized with a pestle, and then the water of the same temperature was emulsified in the lipophilic phase.The cream was stirred continuously, and the active ingredient (2%) was dispersed thoroughly in it.The preservative, propylene glycol [13], was added to the preparation when it cooled to below 30 • C, and stirring was continued until the cream was cooled further to room temperature.The preparation was stored in a refrigerator (2-8 • C) until use.

Viscosity and Particle Size Distribution of Caffeine Cream
Rheological measurements were performed by the Kinexus Pro Rheometer (Malvern Instruments Ltd., Malvern, UK), registering the data with rSpace for Kinexus Pro 1.3 software.Semisolid samples were characterized using a cone and plate geometry (42.5 mm diameter with a 5 • cone) where the gap for sample placement was 0.03 mm.Oscillatory rotational measurements of formulations were determined at 25 • C and controlled with an accuracy of ±0.1 • C by the Peltier system of the instrument.In all measurements (3 parallels) a cylindrical cover made of stainless steel was placed over the samples to create a closed, saturated volume around the sample and prevent evaporation of the samples.
The particle size distribution of the 2% caffeine cream was measured with the Malvern Mastersizer 2000 (Malvern Instruments Ltd., Malvern, UK) laser diffraction device.Then, 1 g of sample was dispersed in 10 mL of water and, after 1 min of Vortex shaking, 1 mL was taken and added to the sample dispenser unit containing 100 mL of purified water.The dispenser unit was applied at 1500 rpm.The multi-channel detector system measured the scattering of red laser light with a wavelength of 633 nm, and the evaluation software provided the results according to the Fraunhofer method after measuring three parallel samples.Each sample was measured in 3 replicates.The average curves were calculated from nine measurements.

Diffusion Platforms
For the assessment of drug permeability, four different barrier models were included.( 1) Polyester (PET) membrane (it4ip, Louvain-la-Neuve, Belgium), which was used since it is a common base of cell culture inserts.(2) Cellulose-acetate membrane with a pore size of 0.45 µm (Sartorius AG, Göttingen, Germany), which was chosen as a control since it is widely used in in vitro release studies [15].( 3) For the third model, we used abdominal rat skin after shaving and depilation followed by mechanical sensitization.Male Wistar rats (ToxiCoop Zrt., Budapest, Hungary, weighing 572-615 g) were housed in an animal room (12 h light/dark periods, 22 ± 3 °C temperature, 50 ± 20% humidity) and fed with commercial laboratory chow and tap water ad libitum.All procedures in this study conformed to the guidelines of the Association for Assessment and Accreditation of Laboratory Animal Care International's expectations for animal use and licensed by the Directorate for the Safety of the Food Chain and Animal Health, Budapest, and Pest County Agricultural Administrative Authority, Hungary (PE/EA/672-6/2016, 8 April 2016).The animals were anesthetized with chloral hydrate (450 mg/kg, i.p.).The abdominal surface was shaved with an electric shaver and depilated with a depilator cream (X-epil, Aveola Kft., Budapest, Hungary).The skin was washed, wiped dry, and mechanically sensitized 10 times by tape stripping with leucoplast (BSN Medical GmbH, Hamburg, Germany) to remove the upper layers of dead keratinocytes.This technique enables detectable concentrations in the dermis and hypodermis.(4) Finally, a blank, acellular hydrogel containing 3% alginate (Merck-Sigma-Aldrich, Darmstadt, Germany) and 1% gelatine (Merck-Sigma-Aldrich, Darmstadt, Germany) was used, which is an important bioink scaffold in tissue engineering.The thicknesses of the various diffusion platforms were the following: PET, CA, rat skin, and alginate scaffold-0.012mm, 0.2 mm, 0.59 mm, and 1.77 mm, measured with a slide caliper (Berger, 020701-0001).The ultrastructure of the synthetic membranes, rat skin, and alginate scaffold were analyzed by scanning electron microscopy (Figure 2).

Diffusion Platforms
For the assessment of drug permeability, four different barrier models were included.(1) Polyester (PET) membrane (it4ip, Louvain-la-Neuve, Belgium), which was used since it is a common base of cell culture inserts.(2) Cellulose-acetate membrane with a pore size of 0.45 µm (Sartorius AG, Göttingen, Germany), which was chosen as a control since it is widely used in in vitro release studies [15].(3) For the third model, we used abdominal rat skin after shaving and depilation followed by mechanical sensitization.Male Wistar rats (ToxiCoop Zrt., Budapest, Hungary, weighing 572-615 g) were housed in an animal room (12 h light/dark periods, 22 ± 3 • C temperature, 50 ± 20% humidity) and fed with commercial laboratory chow and tap water ad libitum.All procedures in this study conformed to the guidelines of the Association for Assessment and Accreditation of Laboratory Animal Care International's expectations for animal use and licensed by the Directorate for the Safety of the Food Chain and Animal Health, Budapest, and Pest County Agricultural Administrative Authority, Hungary (PE/EA/672-6/2016, 8 April 2016).The animals were anesthetized with chloral hydrate (450 mg/kg, i.p.).The abdominal surface was shaved with an electric shaver and depilated with a depilator cream (X-epil, Aveola Kft., Budapest, Hungary).The skin was washed, wiped dry, and mechanically sensitized 10 times by tape stripping with leucoplast (BSN Medical GmbH, Hamburg, Germany) to remove the upper layers of dead keratinocytes.This technique enables detectable concentrations in the dermis and hypodermis.(4) Finally, a blank, acellular hydrogel containing 3% alginate (Merck-Sigma-Aldrich, Darmstadt, Germany) and 1% gelatine (Merck-Sigma-Aldrich, Darmstadt, Germany) was used, which is an important bioink scaffold in tissue engineering.The thicknesses of the various diffusion platforms were the following: PET, CA, rat skin, and alginate scaffold-0.012mm, 0.2 mm, 0.59 mm, and 1.77 mm, measured with a slide caliper (Berger, 020701-0001).The ultrastructure of the synthetic membranes, rat skin, and alginate scaffold were analyzed by scanning electron microscopy (Figure 2).The polymeric structure and porosity are different in cellulose acetate and polyester membranes.The pore size was similar in the two membranes used in this study, but the surface and internal structure were different, as was documented by scanning electron microscopic images.Cellulose acetate (CA) membranes, due to their excellent film-forming property, high chemical and mechanical stability, high hydrophilicity, eco-friendliness, and cost-effectiveness, are extensively used in water and wastewater treatment, gas separation, energy generation, and pharmaceutical and cell culture studies [16].
Polyester (PET) membranes are resistant to most organic solvents, amides, and halogenated hydrocarbons.This broad chemical compatibility makes them suitable for the detection of particles in many corrosive fluids.Also, PE membranes are produced from a pure polymeric film and give exceptional chemical cleanliness.They are free of contaminants and have a low tare weight, minimum water adsorption, and very low levels of nonspecific protein binding.PET membranes are true-pore-size microporous membranes featuring a sharp cut-off, with a closely controlled pore size distribution [17].

Bioanalysis
The caffeine concentrations of the perfused physiological fluids were determined with a UV-VIS spectrophotometer (NanoDrop™ 2000, Thermo Scientific, Budapest, Hungary) after each sample collection.The absorption maximum of caffeine was detected at 272 nm.

Statistical Methods
All data were analyzed and plotted using OriginPro 2015 Software (Northampton, MA, USA).The comparison of experimental groups was carried out using one-way ANOVA followed by the Tukey post hoc test.p values < 0.05 were considered statistically significant.The polymeric structure and porosity are different in cellulose acetate and polyester membranes.The pore size was similar in the two membranes used in this study, but the surface and internal structure were different, as was documented by scanning electron microscopic images.Cellulose acetate (CA) membranes, due to their excellent film-forming property, high chemical and mechanical stability, high hydrophilicity, eco-friendliness, and cost-effectiveness, are extensively used in water and wastewater treatment, gas separation, energy generation, and pharmaceutical and cell culture studies [16].
Polyester (PET) membranes are resistant to most organic solvents, amides, and halogenated hydrocarbons.This broad chemical compatibility makes them suitable for the detection of particles in many corrosive fluids.Also, PE membranes are produced from a pure polymeric film and give exceptional chemical cleanliness.They are free of contaminants and have a low tare weight, minimum water adsorption, and very low levels of nonspecific protein binding.PET membranes are true-pore-size microporous membranes featuring a sharp cut-off, with a closely controlled pore size distribution [17].

Bioanalysis
The caffeine concentrations of the perfused physiological fluids were determined with a UV-VIS spectrophotometer (NanoDrop™ 2000, Thermo Scientific, Budapest, Hungary) after each sample collection.The absorption maximum of caffeine was detected at 272 nm.

Statistical Methods
All data were analyzed and plotted using OriginPro 2015 Software (Northampton, MA, USA).The comparison of experimental groups was carried out using one-way ANOVA followed by the Tukey post hoc test.p values < 0.05 were considered statistically significant.

Computational Fluid Dynamics in Single-Channel Microfluidic Diffusion Chamber (sMDC)
The diffusion behavior of caffeine through permeable membranes or skins within a single-channel microfluidic device was examined by in silico modeling.This study focused on varying flow rates to observe their impact on the accumulation of caffeine and analyze the resultant wall shear stress on the membrane and velocity profiles within the channel.
A single channel of a microfluidic device with one inlet for peripheral perfusion fluid, one diffusion surface inlet for caffeine, and a common outlet were considered, as shown in Figure 3. PPF entered the channel through the inlet at a constant flow rate.Three different flow rates were considered and four different permeable membranes-PET, cellulose acetate, rat skin, and alginate scaffold-were used.They were kept in the section in the middle of the channel.At the top of this membrane, a 2% caffeine cream was applied, which diffused through the permeable membrane into the flow of the PPF.The caffeine-PPF mixture left the domain through the outlet.Figure 3 below shows the geometry of the sMDC.The diffusion behavior of caffeine through permeable membranes or skins within a single-channel microfluidic device was examined by in silico modeling.This study focused on varying flow rates to observe their impact on the accumulation of caffeine and analyze the resultant wall shear stress on the membrane and velocity profiles within the channel.
A single channel of a microfluidic device with one inlet for peripheral perfusion fluid, one diffusion surface inlet for caffeine, and a common outlet were considered, as shown in Figure 3. PPF entered the channel through the inlet at a constant flow rate.Three different flow rates were considered and four different permeable membranes-PET, cellulose acetate, rat skin, and alginate scaffold-were used.They were kept in the section in the middle of the channel.At the top of this membrane, a 2% caffeine cream was applied, which diffused through the permeable membrane into the flow of the PPF.The caffeine-PPF mixture left the domain through the outlet.Figure 3 below shows the geometry of the sMDC.A finer version of this mesh was used for simulations.

Computational Fluid Dynamics in Single-Channel Microfluidic Diffusion Chamber (sMDC)
Model used for modeling fluid flow (cream penetration) through the porous layers: The species transport model was used to model the diffusion process with the inletdiffusion option and the energy equation turned on.
Governing equations: For the perfusion fluid: Continuity Equation [18]: Momentum Equation: Species Transport Equation [19]: where Y i is the mass fraction of the i th species and → J i is the diffusion flux of the i th species, which is given by where D i is the mass diffusion constant of the i th species.
For the flow through porous layers [20]: Continuity Equation: Momentum Equation: Species Transport Equation [19]: In the above Equations ( 1)-( 5), ∇ is the gradient operator, ρ is the density of the fluid perfusion fluid, and u = u x i+u y j,+ u z k, is the velocity vector, with u x , u y , and u z as its respective components in the x (horizontal), y (vertical), and z directions.Du Dt (= ∂u ∂t + u•∇u) is the material derivative of the velocity vector u, representing the total (local plus convective) change in velocity experienced by a fluid particle.Pressure, time, and body force (e.g., gravity) follow their standard symbols as p, t, and f, respectively.The perfusion fluid viscosity is represented by µ.The properties of the membranes are designated as porosity by ϵ, permeability by κ, and Forchheimer inertial factor by F.

Numerical Method
A computational fluid dynamics tool was used-Ansys Fluent [21]-for simulating the flow of PPF and the diffusion of caffeine creams into PET, CA, alginate, and rat skin.It was assumed that PPF is a Newtonian fluid with constant viscosity.Since the membranes were porous in nature, the porous media model was used for modeling flow through these barriers.Accordingly, the governing equations were the Navier-Stokes equation, the concentration equation in the channel, the Darcy-Brinkman-Forchheimer model [20], and the modified concentration equation for the four membranes.

Parameters Used
The simulation study investigated flow rates of PPF at 4, 40, and 100 µL/min, using an estimated viscosity of 0.0009012 Pa•s for PPF [22,23].Caffeine concentration was maintained at 2%.Four different membranes-PET, CA, rat skin, and alginate scaffold-were tested at these flow rates, with thicknesses of 0.012 mm, 0.2 mm, 0.59 mm, and 1.77 mm, respectively.This study explored various Darcy numbers (10 −4 , 10 −6 , 10 −8 ), with a Darcy number of 10 −8 or higher preventing PPF penetration into the membranes, allowing only diffusion and leading to a converged solution.Therefore, Darcy numbers of 10 −8 and higher were utilized.
It was assumed that the PPF was a Newtonian fluid with constant viscosity, and the viscosity of the mixture of caffeine and PPF was allowed to be varied according to the mass-weighted viscosity law [24].However, due to a very low mass fraction (in the order of 10 −5 or 10 −6 ) of caffeine in the flow during steady state, the viscosity of the mixture remained nearly constant near the value of PPF.Thus, the viscosity of caffeine did not have a significant effect on the simulations.

Simulations
The simulation process involved several steps.Initially, structured mesh files corresponding to different membranes were created using Ansys Workbench.These mesh files were then imported into Fluent for flow analysis.Variation in flow rate was achieved by adjusting the velocity at the inlet of the PPF.
Due to a lack of data about the accurate mass diffusion coefficient of caffeine into the perfusion fluid, in order to analyze the diffusive behavior of caffeine, its diffusive velocity was iteratively determined at the upper surface of the porous membrane to obtain the diffusion rate that agreed with the experimental results.Unsteady simulations were executed with a timestep of 3 s, spanning a total duration of 18,000 s (equivalent to 5 h).
At the conclusion of each timestep, the caffeine mass flux was estimated by the following formula: . m m, t is the mass flow rate of the mixture at timestep t, Y c, t is the mass fraction of caffeine, and A is the area through which diffusion occurs.Subsequently, this value was numerically integrated over time to compute the accumulated mass per diffusion area.

Viscosity and Particle Size Distribution
The particle size distribution study showed a homogenous composition of the 2% caffeine cream, which is presented in Figure 4A,B.In total, 90% of the particles appeared to be 140 µm in size and only one peak could be seen in the frequency curve.The viscosity of the cream was determined at different shear rates from 0.1 to 100/s (Figure 4C), and it showed an almost linear inverse relationship with the increasing shear rate (from 800 to 2 Pa• s) due to the shear-thinning behavior of the cream.The release and diffusion of caffeine from the cream formulation and across the diffusion platforms (PET and CA membranes, excised rat skin, and alginate scaffold) were investigated at three different flow rates (4, 40, and 100 µL/min) for 5 h.Three diffusion chambers were compared (sMDC: single-channel microfluidic diffusion chamber, The release and diffusion of caffeine from the cream formulation and across the diffusion platforms (PET and CA membranes, excised rat skin, and alginate scaffold) were investigated at three different flow rates (4, 40, and 100 µL/min) for 5 h.Three diffusion chambers were compared (sMDC: single-channel microfluidic diffusion chamber, mMDC: multi-channel microfluidic diffusion chamber, and LiveBox2).The cumulative mass-time profiles are displayed in Figures 5-8.
Figure 5 displays the caffeine permeation profile across the polyester (PET) membrane in different setups: sMDC, mMDC, and LiveBox2.Across all setups, permeation at 40 µL/min was the highest (AUC was 1.7-, 1.8-, and 12.9-fold compared to the 4 µL/min flow rate in sMDC, mMDC, and Livebox2, respectively).However, at a flow rate of 100 µL/min, permeation decreased compared to the flow rate of 40 µL/min.In Figure 6, the permeation profile of caffeine through the cellulose-acetate membrane is depicted.The total penetrated caffeine amount was markedly higher in sMDC and LiveBox2 at a flow rate of 40 µL/min than at 4 µL/min (AUC was 5.9-fold and 13-fold, respectively).This increase was attributed to the low caffeine concentration in the receptor chamber at higher flow rates, leading to an elevated concentration gradient that facilitated drug diffusion.Despite a greater concentration difference at a 100 µL/min flow rate, the cumulative mass increased only in the LiveBox2 system, although the difference was not statistically significant.The absorption parameters of caffeine through rat skin, as shown in Figure 7, were lower compared to the previous barrier models (e.g., PET and CA membranes).Here, also, a significant difference was observed between the first two flow rates (4 and 40 µL/min) in all setups, but the total penetrated caffeine decreased at 100 µL/min.Figure 8 depicts the results obtained with the alginate scaffold.The highest caffeine levels were recorded at 40 µL/min (AUC was 1.4, 2, and 16 times greater than at the lowest fluid velocity), but they decreased at the highest flow rate.Figure 5 displays the caffeine permeation profile across the polyester (PET) membrane in different setups: sMDC, mMDC, and LiveBox2.Across all setups, permeation at 40 µL/min was the highest (AUC was 1.7-, 1.8-, and 12.9-fold compared to the 4 µL/min flow rate in sMDC, mMDC, and Livebox2, respectively).However, at a flow rate of 100 µL/min, permeation decreased compared to the flow rate of 40 µL/min.In Figure 6, the permeation profile of caffeine through the cellulose-acetate membrane is depicted.The total penetrated caffeine amount was markedly higher in sMDC and LiveBox2 at a flow rate of 40 µL/min than at 4 µL/min (AUC was 5.9-fold and 13-fold, respectively).This    chamber size was the smallest, the differences in the penetration profiles were less prominent.Significant variations were observed among barrier models, too.In sMDC, the highest penetration rate occurred in the CA membrane, followed by PET, alginate, and excised rat skin.In LiveBox2, the highest permeability was in the alginate scaffold, followed by the CA membrane, polyester membrane, and rat skin.Conversely, in mMDC, all barrier models exhibited similar ranges of drug penetration.Based on the penetration results, the molecular diffusion greatly depended on the geometric design of the microfluidic device.In the case of the single-channel device, for the PET and CA membranes and also the rat skins, the 40 µL/min flow rate produced the highest diffusion, while, in LiveBox2, the highest flow rate (100 µL/min) led to the most effective diffusion for the CA membranes and rat skin.Thus, similarly to blood circulation, in the LiveBox2 system, the channel diameter was larger than in sMDC and, therefore, the optimal flow rate was faster than in the case of the microchannel systems.
All in all, the medium flow rate (40 µL/min) seemed to be the optimum perfusion velocity in the majority of the diffusion platforms tested in these microfluidic setups.Both the slower and the faster flow rates showed worse diffusion effectivity for caffeine through the barriers.The best barrier was the natural animal skin, which displayed the lowest The LiveBox2 system possessed the highest receptor chamber volume, leading to the largest disparities in caffeine permeation between the first two flow rates.In contrast, the sMDC was designed with a considerably smaller volume, resulting in less pronounced differences between results at 4 µL/min and 40 µL/min.Similarly, in mMDC, where the chamber size was the smallest, the differences in the penetration profiles were less prominent.
Significant variations were observed among barrier models, too.In sMDC, the highest penetration rate occurred in the CA membrane, followed by PET, alginate, and excised rat skin.In LiveBox2, the highest permeability was in the alginate scaffold, followed by the CA membrane, polyester membrane, and rat skin.Conversely, in mMDC, all barrier models exhibited similar ranges of drug penetration.
Based on the penetration results, the molecular diffusion greatly depended on the geometric design of the microfluidic device.In the case of the single-channel device, for the PET and CA membranes and also the rat skins, the 40 µL/min flow rate produced the highest diffusion, while, in LiveBox2, the highest flow rate (100 µL/min) led to the most effective diffusion for the CA membranes and rat skin.Thus, similarly to blood circulation, in the LiveBox2 system, the channel diameter was larger than in sMDC and, therefore, the optimal flow rate was faster than in the case of the microchannel systems.
All in all, the medium flow rate (40 µL/min) seemed to be the optimum perfusion velocity in the majority of the diffusion platforms tested in these microfluidic setups.Both the slower and the faster flow rates showed worse diffusion effectivity for caffeine through the barriers.The best barrier was the natural animal skin, which displayed the lowest permeability for the test compound due to the presence of tight junction proteins between the keratinocytes in the epidermal stratum granulosum layer, the presence of active transporter proteins in the skin cells, and the complex extracellular matrix.The CA membrane and alginate produced the highest diffusion (permeability), which peaked with sMDC for the CA membrane and with LiveBox2 for the alginate.In sMDC and mMDC, the microfluidic channels were narrow and the artificial membranes showed hydrophilic characteristics.On the other hand, in the skin, the penetration was more moderate due to the complexity of the natural epidermal/dermal barrier, including the lipid matrix in the stratum corneum, which presented a lipophilic surface, partially active transporters, and tight junctions proteins in the granular epidermis.Caffeine MW, LogP, particle size, viscosity, and the geometry of the device altogether resulted in the best balance for the diffusion at a 40 µL/min flow rate in these devices in almost all cases investigated.

Computational Fluid Dynamics (CFD)
The CFD simulations were only performed for the sMDC microfluidic device.The results of the simulation included figures depicting the cumulative mass versus time (Figure S1), non-dimensional velocity profiles (Figure S2) at various cross-sections along the duct's length, caffeine concentration contours (Figure 12), and skin wall shear contours (Figure 14).The cumulative mass plots exhibited nearly linear trends after a certain duration, indicating a nearly constant diffusion rate at a steady state (Figure S1).Lower flow rates required more time to reach a steady state, whereas higher flow rates achieved a steady state within approximately 10 min.
The contours of shear stress at a flow rate of 40 µL/min are presented in Figure 9A-D for PET, CA, rat skin, and alginate scaffold.In all the figures, some amount of stress concentration near the beginning and end of the membrane can be seen.This can be explained by the sudden change in the cross-sectional area at these regions, which caused the flow profile to change by a small amount.This change resulted in stress being concentrated in these regions.It can also be observed that only the base of the membrane suffered shear stress (Figure 9B-D), which was because that was the only surface exposed to the flow of PPF owing to the low permeability of the membranes [25].However, due to the very low thickness of the PET membrane, the shear can be seen everywhere in the membrane in Figure 9A.permeability for the test compound due to the presence of tight junction proteins between the keratinocytes in the epidermal stratum granulosum layer, the presence of active transporter proteins in the skin cells, and the complex extracellular matrix.The CA membrane and alginate produced the highest diffusion (permeability), which peaked with sMDC for the CA membrane and with LiveBox2 for the alginate.In sMDC and mMDC, the microfluidic channels were narrow and the artificial membranes showed hydrophilic characteristics.On the other hand, in the skin, the penetration was more moderate due to the complexity of the natural epidermal/dermal barrier, including the lipid matrix in the stratum corneum, which presented a lipophilic surface, partially active transporters, and tight junctions proteins in the granular epidermis.Caffeine MW, LogP, particle size, viscosity, and the geometry of the device altogether resulted in the best balance for the diffusion at a 40 µL/min flow rate in these devices in almost all cases investigated.

Computational Fluid Dynamics (CFD)
The CFD simulations were only performed for the sMDC microfluidic device.
The results of the simulation included figures depicting the cumulative mass versus time (Figure S1), non-dimensional velocity profiles (Figure S2) at various cross-sections along the duct's length, caffeine concentration contours (Figure 12), and skin wall shear contours (Figure 14).The cumulative mass plots exhibited nearly linear trends after a certain duration, indicating a nearly constant diffusion rate at a steady state (Figure S1).Lower flow rates required more time to reach a steady state, whereas higher flow rates achieved a steady state within approximately 10 min.
The contours of shear stress at a flow rate of 40 µL/min are presented in Figure 9A-D for PET, CA, rat skin, and alginate scaffold.In all the figures, some amount of stress concentration near the beginning and end of the membrane can be seen.This can be explained by the sudden change in the cross-sectional area at these regions, which caused the flow profile to change by a small amount.This change resulted in stress being concentrated in these regions.It can also be observed that only the base of the membrane suffered shear stress (Figure 9B-D), which was because that was the only surface exposed to the flow of PPF owing to the low permeability of the membranes [25].However, due to the very low thickness of the PET membrane, the shear can be seen everywhere in the membrane in Figure 9A.

Shear Stress Profiles of Membranes
Variations of shear stress along the bottom surface of the PET, CA, rat skin, and alginate scaffold are shown in Figure 10.It can be seen that the maximum shear stress at different flow rates had different values across the membranes.PET had higher values of shear stress than the other membranes, as seen in Figure 10A; this was due to its low thickness [25], offering less resistance to shear stress.It can also be observed in Figure  [22,23]).A flow rate of 40 µL/min was selected as the most effective flow for caffeine diffusion.

Shear Stress Profiles of Membranes
Variations of shear stress along the bottom surface of the PET, CA, rat skin, and alginate scaffold are shown in Figure 10.It can be seen that the maximum shear stress at different flow rates had different values across the membranes.PET had higher values of shear stress than the other membranes, as seen in Figure 10A; this was due to its low thickness [25], offering less resistance to shear stress.It can also be observed in Figure 10A,D that the alginate had the highest shear stress at a flow rate of 40 µL/min, while PET had the highest shear stress at a flow rate of 100 µL/min.Cellulose acetate experienced the least amount of shear stress at all flow rates, which is clearly visible in Figure 10B.[22,23]), (C skin (h = 0.59 mm, Da = 10 −8 , porosity = 0.03 (assumed)), and (D) alginate scaffold (h = 1.77 mm = 10 −10 , porosity = 0.4 [22,23]).A flow rate of 40 µL/min was selected as the most effective flow caffeine diffusion.

Shear Stress Profiles of Membranes
Variations of shear stress along the bottom surface of the PET, CA, rat skin, and a nate scaffold are shown in Figure 10.It can be seen that the maximum shear stress at ferent flow rates had different values across the membranes.PET had higher value shear stress than the other membranes, as seen in Figure 10A; this was due to its l thickness [25], offering less resistance to shear stress.It can also be observed in Fig 10A,D that the alginate had the highest shear stress at a flow rate of 40 µL/min, while P had the highest shear stress at a flow rate of 100 µL/min.Cellulose acetate experienced least amount of shear stress at all flow rates, which is clearly visible in Figure 10B.

Velocity Profiles
The velocity profiles below the various membranes were also calculated in the microchannel of sMDC and the results are displayed in Figure S2.It shows the non-dimensional velocity profiles along the height of the channel at different locations below the membranes.The origin was taken to be in the middle of the channel at the bottom wall.The velocity profiles were then plotted at the symmetrically opposite points: x = ±0.02m indicated points nearly midway between the inlet and membrane and the membrane and outlet; x = ±0.00445m indicated the beginning and end of the membrane; x = 0 indicated the point exactly in the middle of the chamber.The symmetrically opposite points had the same velocity profiles and, thus, they overlapped.The little difference between the profiles at x = 0, ±0.00445 and ±0.02 could be explained by the existence of a membrane that was not a perfect wall and caused the flow to deviate a little from its fully developed profiles.At the same time, the membranes had very low permeability; thus, different membranes did not cause an observable change in the velocity profiles, which can be seen in the very similar profiles for each membrane.

Velocity Contours
Figure 11A-C show the contours of velocity for different flow rates for all four membranes.We can see that no flow penetrated any membrane at any flow rate.Thus, the flow profile did not change significantly throughout the channel and the contours in the channel were nearly the same for all the membranes at a given flow rate.This can be explained by the fact that the membranes were almost impervious to flow due to their low permeability; thus, flow in the channel behaved similarly to Plane-Poiseuille flow [26].
same velocity profiles and, thus, they overlapped.The little difference between the profi at x = 0, ±0.00445 and ±0.02 could be explained by the existence of a membrane that w not a perfect wall and caused the flow to deviate a little from its fully developed profi At the same time, the membranes had very low permeability; thus, different membra did not cause an observable change in the velocity profiles, which can be seen in the v similar profiles for each membrane.

Velocity Contours
Figure 11A-C show the contours of velocity for different flow rates for all four me branes.We can see that no flow penetrated any membrane at any flow rate.Thus, the fl profile did not change significantly throughout the channel and the contours in the ch nel were nearly the same for all the membranes at a given flow rate.This can be explain by the fact that the membranes were almost impervious to flow due to their low perm bility; thus, flow in the channel behaved similarly to Plane-Poiseuille flow [26].

Caffeine Progression
In Figure 12A-C, the time evolution of drug diffusion through rat skin into the perfusion fluid (PPF) is shown for flow rates of 4 µL/min, 40 µL/min, and 100 µL/min, respectively.At a flow rate of 4 µL/min, the system approached a steady state gradually.For the first 1800 s, almost no mass passed through the outlet, as indicated by the horizontal line in the cumulative mass graph (Figure 7A) during the initial 30 min.After this period, the system reached a steady state, resulting in a nearly constant mass flux and a linear increase in caffeine accumulation.The concentration profiles were also thicker, as seen in the contours.
At a higher flow rate of 40 µL/min, drug diffusion followed a similar pattern but reached a steady state much faster, in just 200 s (Figure 12B).A notable thick concentration layer formed on the top wall of the channel due to the higher flow velocities.At the highest flow rate of 100 µL/min, a steady state was achieved in a similar timeframe.However, the concentration remained very low, with barely any mass passing through the outlet, as shown in Figure 12C.
These results illustrate how flow rate significantly influenced the time to reach a steady state and the concentration profiles in the system.Lower flow rates took longer to reach a steady state but resulted in higher concentrations, while higher flow rates achieved a steady state quickly but with much lower concentrations.The higher concentrations may have been due to the long exposure time for lower flow rates, which allowed more mass of caffeine to accumulate in the flow.
These results illustrate how flow rate significantly influenced the time to rea steady state and the concentration profiles in the system.Lower flow rates took long reach a steady state but resulted in higher concentrations, while higher flow rates achi a steady state quickly but with much lower concentrations.The higher concentrations have been due to the long exposure time for lower flow rates, which allowed more m of caffeine to accumulate in the flow.

Three-Dimensional Simulations
Two-dimensional simulation assumed that the channel maintained a constant c section throughout the flow.However, in reality, the cross-section changed in the y a directions near the diffusion area.The effect of the change in the cross-section in th direction is discussed in Sections 3.3 and 3.5 and visible in Figures 9 and S2.To exam the impact of the cross-sectional change in the z-direction, 3D results for alginate at a rate of 4 µL/min are presented.The streamlines shown in Figure 13 indicate that there no significant flow penetration into the membrane.
Figure 14 reveals an important result that was not visible in the 2D simulation.cross-section of the duct changed in the z-direction in the middle of the channel, w diffusion occurred.This change caused variations in velocity in this region that wer evident in the 2D velocity contours (Figure 11).Consequently, there was a gradie velocity here that could not be visualized in the 2D simulation, leading to stress con trations (in addition to the ones observed in Figure 9) at the beginning and end o membrane (Figure 14).This increase in stress is not apparent in the shear profiles sh earlier (Figure 9).

Three-Dimensional Simulations
Two-dimensional simulation assumed that the channel maintained a constant crosssection throughout the flow.However, in reality, the cross-section changed in the y and z-directions near the diffusion area.The effect of the change in the cross-section in the y-direction is discussed in Sections 3.3 and 3.5 and visible in Figures 9 and S2.To examine the impact of the cross-sectional change in the z-direction, 3D results for alginate at a flow rate of 4 µL/min are presented.The streamlines shown in Figure 13 indicate that there was no significant flow penetration into the membrane.These findings highlight the limitations of 2D simulations in capturing the full complexity of the system.The 3D simulations provide a more accurate representation by accounting for cross-sectional changes and their effects on flow dynamics.The future scope for this work could involve exploring 3D modeling of these channels for a better and more in-depth understanding of the flow pattern near the diffusion area.

Discussion and Conclusions
A state-of-the-art method of computational fluid dynamics has been successfully applied for modeling fluid flow in microfluidic systems with different complexities.Bakuova and co-workers investigated and compared the flow behavior and filling characteristics of two microfluidic liver-on-a-chip devices using CFD analysis and experimental cell cul- Figure 14 reveals an important result that was not visible in the 2D simulation.The cross-section of the duct changed in the z-direction in the middle of the channel, where diffusion occurred.This change caused variations in velocity in this region that were not evident in the 2D velocity contours (Figure 11).Consequently, there was a gradient of velocity here that could not be visualized in the 2D simulation, leading to stress concentrations (in addition to the ones observed in Figure 9) at the beginning and end of the membrane (Figure 14).This increase in stress is not apparent in the shear profiles shown earlier (Figure 9).These findings highlight the limitations of 2D simulations in capturing the full complexity of the system.The 3D simulations provide a more accurate representation by accounting for cross-sectional changes and their effects on flow dynamics.The future scope for this work could involve exploring 3D modeling of these channels for a better and more in-depth understanding of the flow pattern near the diffusion area.

Discussion and Conclusions
A state-of-the-art method of computational fluid dynamics has been successfully applied for modeling fluid flow in microfluidic systems with different complexities.Bakuova and co-workers investigated and compared the flow behavior and filling characteristics of two microfluidic liver-on-a-chip devices using CFD analysis and experimental cell culture growth based on the Huh7 cell line [27].Biagini and co-workers developed a millifluidic device to generate fully controlled shear stress profiles and quantitatively probe their influence on tissue or bacterial models, overcoming the limitations of previously reported similar devices [28].Amini and his group studied the flow behavior of a liquid-liquid (chloroform and water) extraction process in a serpentine microchannel [29].The simulation was performed using a 3D model and the were found to be consistent with experimental data.Cardona and her research group used numerical modeling for physical cell trapping in microfluidic chips.The validity of the model was assessed with experimental data [30].In a recent paper, Garud and co-workers applied a hybrid approach of computational fluid dynamics and Taguchi analysis to evaluate the influence of various geometrical and operating factors on shear stress in a microfluidic cerebrovascular channel [31].Drug delivery across the blood-brain barrier was considered in the presence of

Discussion and Conclusions
A state-of-the-art method of computational fluid dynamics has been successfully applied for modeling fluid flow in microfluidic systems with different complexities.Bakuova and co-workers investigated and compared the flow behavior and filling characteristics of two microfluidic liver-on-a-chip devices using CFD analysis and experimental cell culture growth based on the Huh7 cell line [27].Biagini and co-workers developed a millifluidic device to generate fully controlled shear stress profiles and quantitatively probe their influence on tissue or bacterial models, overcoming the limitations of previously reported similar devices [28].Amini and his group studied the flow behavior of a liquid-liquid (chloroform and water) extraction process in a serpentine microchannel [29].The simulation was performed using a 3D model and the results were found to be consistent with experimental data.Cardona and her research group used numerical modeling for physical cell trapping in microfluidic chips.The validity of the model was assessed with experimental data [30].In a recent paper, Garud and co-workers applied a hybrid approach of computational fluid dynamics and Taguchi analysis to evaluate the influence of various geometrical and operating factors on shear stress in a microfluidic cerebrovascular channel [31].Drug delivery across the blood-brain barrier was considered in the presence of shear stress at the blood-brain interface.Another vascular chip model was analyzed by Wang's group [32].CFD models were employed to reveal the hydrodynamics in a tissue-engineered blood vessel on a chip for cardiovascular disease examinations.Takken and Wille described an accelerated CFD simulation method of microfluidic devices by exploiting higher levels of abstraction [33].They showed case studies that confirm that the proposed method accelerates CFD simulations by multiple factors (often by several orders of magnitude) while maintaining the fidelity of the simulations.
In our report, a novel application field of CFD in the development and characterization of microfluidic skin-on-a-chip devices was described.The core idea of the current study was to analyze and optimize the fluid dynamics in three existing microfluidic diffusion chambers for transdermal drug delivery studies.Due to the limitation of the analytical detection techniques of active ingredients in the perfusion fluid in such microscale systems, the aim was to find the highest extracting flow rate and obtain measurable drug concentrations in the collected perfusate samples.This study included (1) an experimental part, where three perfusion flow rates were studied using four diffusion platforms integrated into different microfluidic diffusion chambers, and (2) a 2D and 3D simulation part, in which we wanted to find the theoretical explanation for our experimental findings.
The experimental part revealed that the highest drug diffusion could be achieved at the medium flow rate (40 µL/min) in all three microfluidic setups in almost all experimental arrangements, and the barrier function of the natural skin was the strongest, while the permeability turned out to be greatest in the CA membrane in sMDC and the alginate scaffold in the LiveBox2 system.In terms of extracting the absorbed active ingredient, the mMDC system was the least efficient.The smallest cumulative amounts were measured here.In terms of ease of application, the sMDC system was the best to handle, while the LiveBox2 chamber seemed to be the most advantageous from a practical point of view in terms of washing and cleanability (including sterilization).The experimental findings show the direction for further developments of microfluidic systems for drug penetration testing.Based on the results, sMDC and mMDC should be used at lower flow rates, especially in the case of hydrophilic drugs, while LiveBox2 works better at higher speeds.This study also showed preliminary data about alginate as a scaffold bioink for the production of 3D-bioprinted skin tissue for transdermal drug absorption testing.Simulations should be continued for the LiveBox 2 system with bioprinted skin.
The in silico modeling showed that the shear rate was the highest at the highest flow rate (100 µL/min) in each device, and, at this velocity, the penetration through the excised rat skin was very low.Although the concentration gradient between the barrier surface and the flowing fluid was proportional to the perfusion flow rate, the degree of diffusion could not be enhanced further due to the inhibitory shear forces.The simulation nicely visualized the shear stress areas of the channels and the caffeine progression with the peripheral perfusion fluid through the microchannel up to the outlet on a longitudinal scale.
Computational fluid dynamics simulation played a crucial role in analyzing a singlechannel microfluidic system.It accurately predicted fluid flow patterns within the microchannel, essential for understanding fluid behavior.CFD also provided detailed insights into shear stress distribution along membrane surfaces, highlighting areas of high stress concentration that could impact system performance.By visualizing concentration gradients, such as caffeine diffusion, how different flow rates influenced the process was revealed.The simulations were validated against experimental data, confirming their reliability.Also, considering three-dimensional effects, such as changes in cross-sectional areas and velocity gradients, offered a more comprehensive understanding compared to two-dimensional analysis.Overall, CFD simulations were indispensable for predicting fluid behavior, validating experimental results, and understanding complex interactions within the microfluidic system.
The consistency between the visualization (CFD simulations) and measured data underscores the reliability and robustness of the microfluidic system.This alignment suggests that the CFD models accurately capture the physical phenomena occurring within the microfluidic device, such as fluid flow patterns, shear stress distributions, and concentration gradients.Specifically, the ability of the simulations to replicate the experimentally observed caffeine diffusion and shear stress patterns validates the precision of the computational models.This validation enhances confidence in the system's design and its predictive capabilities, indicating that the microfluidic platform can reliably simulate real-world conditions for drug delivery studies and other applications.Moreover, the accurate reflection of three-dimensional effects and localized phenomena in the simulations further confirms the system's efficacy in providing detailed and accurate representations of microscale fluid dynamics, thereby supporting its use in designing and optimizing microfluidic devices for various biomedical applications.
The computational investigation into caffeine diffusion through a permeable membrane within a microfluidic device revealed several significant findings.The cumulative mass plots illustrated a near-linear trend after an initial period, indicative of a steady diffusion rate once equilibrium was reached.Notably, lower flow rates necessitated a longer duration to attain a steady state, while higher flow rates achieved equilibrium within approximately 10 min.Stress concentration near the beginning and end of the membrane was observed across all figures, attributed to abrupt changes in the cross-sectional area leading to localized alterations in the flow profile and subsequent stress concentration.Interestingly, the distribution of shear stress along the membrane surfaces varied, with the PET membrane exhibiting higher shear values compared to other membranes due to its low thickness, which offered less resistance to shear.Additionally, the shear distribution varied with flow rate, with alginate showing the highest shear stress among membranes

Figure 3 .Figure 3 .
Figure 3. (A,B) Schematic of the problem, (C) Representation of fluid and porous regions in the computational domain, (D) Coarse mesh (uniform mesh), (E) Fine mesh (finer in the critical regions).A finer version of this mesh was used for simulations.2.7.1.Computational Fluid Dynamics in Single-Channel Microfluidic Diffusion Chamber (sMDC) Model used for modeling fluid flow (cream penetration) through the porous layers: The species transport model was used to model the diffusion process with the inletdiffusion option and the energy equation turned on.Governing equations: For the perfusion fluid: Sci. Pharm.2024, 92, x FOR PEER REVIEW 9 of 22

Figure 14 .
Figure 14.Shear stress below alginate scaffold in sMDC device at Q = 4 µL/min: (A) front view; (B) magnified front view near the alginate membrane; (C) bottom view; (D) shear stress at alginate exposed to the flow of PPF.

Figure 14 .
Figure 14.Shear stress below alginate scaffold in sMDC device at Q = 4 µL/min: (A) front view; (B) magnified front view near the alginate membrane; (C) bottom view; (D) shear stress at alginate exposed to the flow of PPF.

Figure 14 .
Figure 14.Shear stress below alginate scaffold in sMDC device at Q = 4 µL/min: (A) front view; (B) magnified front view near the alginate membrane; (C) bottom view; (D) shear stress at alginate exposed to the flow of PPF.

Table 1 .
The final excipient composition of 2% caffeine creams.

Table 2 .
Technical details of the different diffusion systems.sMDC: single-channel microfluidic diffusion cell, mMDC: multichannel microfluidic diffusion cell.