Computational Study of pH-sensitive Hydrogel-based Microfluidic Flow Controllers

This computational study investigates the sensing and actuating behavior of a pH-sensitive hydrogel-based microfluidic flow controller. This hydrogel-based flow controller has inherent advantage in its unique stimuli-sensitive properties, removing the need for an external power supply. The predicted swelling behavior the hydrogel is validated with steady-state and transient experiments. We then demonstrate how the model is implemented to study the sensing and actuating behavior of hydrogels for different microfluidic flow channel/hydrogel configurations: e.g., for flow in a T-junction with single and multiple hydrogels. In short, the results suggest that the response of the hydrogel-based flow controller is slow. Therefore, two strategies to improve the response rate of the hydrogels are proposed and demonstrated. Finally, we highlight that the model can be extended to include other stimuli-responsive hydrogels such as thermo-, electric-, and glucose-sensitive hydrogels.


Introduction
As a relatively new branch of science and technology, microfluidics, which emerged in the early 1990s [1][2][3], has attracted much attention for its diverse applications, ranging from ink-jet printers and fuel injection [4], over surface processing and biological assay [5], to control system, heat management and display technology [6]. Other fields where microfluidic systems are considered and employed include micromixing [7,8], biology and biochemical analysis [9][10][11][12]. For several of these applications, the ability to manipulate the fluid flow within the microchannels is essential [12]; therefore, considerable effort has been devoted to develop microfluidic flow controllers. The majority of flow controllers in microfluidics systems are miniaturized version of their conventional macroscale counterparts [13], which are generally integrated devices comprising electrical, mechanical and optical elements with individual functions. These conventional microfluidic flow controllers have two major drawbacks: the inherent difficulty in assembling the various components into a single system and the requirement of an external power supply, both of which limit their implementation in numerous applications [14,15].
In contrast to conventional microfluidic flow controllers, stimuli-sensitive hydrogels can be employed without external power supply. Moreover, hydrogels offer significant reduction in the complexity of a microsystem due to their unique stimuli-sensitive ability; that is, stimuli-sensitive hydrogels can sense changes in its environment-temperature, pH, glucose, electric field and pressure-and then swell or shrink correspondingly [16,17]. During swelling, certain hydrogels are able to absorb large amounts of water leading to a large swelling ratio [18][19][20]. Hydrogels could therefore replace the major components in microfluidics flow controllers such as sensors, signal processors, regulators and actuators [21]. Moreover, their high water content and soft consistencies lend them excellent biocompatibility, allowing application of this hydrogel-based system in biomedical and biotechnical fields. Due to their potential for active flow control, numerous experimental designs and studies of hydrogels as flow controller in microfluidic systems have been conducted [15,[21][22][23][24][25][26][27]. In contrast, few studies concerning mathematical modeling and simulation have been reported [28][29][30]; hence, it is of interest to develop mathematical models which can aid in the synthesis and design of hydrogels as microfluidic flow controllers. This paper addresses the sensing and actuating behavior of hydrogels as microfluidic flow controller with the aim to derive and analyze a simple mathematical model for a pH-sensitive hydrogel that can be integrated with the external flow in a microfluidic flow system; e.g., in a T-junction, as illustrated in Figure 1. The mathematical model, which takes into account conservation of momentum, mass and ions for laminar incompressible flow and the sensing/deformation of a pH-sensitive hydrogel, is derived, analyzed and presented in Section 2. Details of the numerical procedure are outlined in Section 3. Calibration and validation with steady-state and transient experiments [31] is then carried out for the deformation of the hydrogel as a function of pH, after which we demonstrate how the model can be employed to study the sensing and actuating behavior of the hydrogel as microfluidic flow controller. The flow configurations considered in this paper are (i) a T-junction with hydrogels in one branch and (ii) a T-junction with hydrogels in each branch; the latter has two hydrogels with opposite behavior: a positive pH-responsive hydrogel, which swells as the pH increases, and a negative pH-responsive hydrogel, which shrinks as the pH increases. We finish with conclusions, in which we highlight how the model can be generalized for other types of stimuli-responsive hydrogels.

Mathematical Formulation
In this section, we derive a mathematical model that incorporates the conservation of momentum, mass and species for a laminar incompressible flow as well as the sensing and deformation of the hydrogel. The hydrogel considered in this paper is pHEMA (polyhydroxylethylmethacrylate), which is a pH-sensitive hydrogel [31]. The hydrogel is embedded as microfluidic flow controller in a microchannel (see Figure 1); as the hydrogel shrinks and swells depending on the pH of the solution in the system (sensing), it affects (controls) the overall flow in the system. The solution is aqueous with protons (H + ), sodium ions (Na + ), hydroxide ions (OH − ), and chloride ions (Cl − ) at ambient temperature.

Governing Equations
For laminar, incompressible flow inside the microchannel, conservation of mass and momentum are given by the Navier-Stokes equation: where ρ We account for ion transfer inside the channel and the hydrogel with the Nernst-Planck equation and electroneutrality condition [32][33][34], which can be expressed as where F is Faraday's constant, R is the universal gas constant, ψ is the electric potential, D k is the diffusive coefficient, z k and c k are valence and concentration of the ion species k (= H + , Na + , OH − and Cl − ) respectively; 0 is the permittivity for vacuum, is the dielectric constant of medium relative to vacuum, z f is the valence of fixed charge and c f is the fixed charge concentration inside the hydrogel. For the hydrogel, conservation of mass is solved for a biphasic mixture comprising the solid and fluid phase, whereas conservation of momentum for the hydrogel is considered in terms of Navier's equation with infinitesimal deformations and generalized Darcy's law for a moving porous medium [35][36][37]; that is where q =ρ (f ) v (r) is the Eulerian relative flow vector of the fluid phase with respect to the polymer is the velocity of fluid relative to the solid phase velocity, v (i) is the intrinsic velocity of phase i (solid and fluid), κ (p) is the permeability of polymer phase, σ is the mixture stress tensor, and 0 is the true density of each phase. The velocity of the solid phase can be defined as the rate of deformation; that is, v (p) = ∂u (p) /∂t, where u is the deformation of the hydrogel.

Constitutive Relations
The mixture stress tensor for the hydrogel is given by [35] ef f is the elastic stress tensor of the polymer phase in the hydrogel. We treat the polymer phase as an isotropic elastic material, whence the elastic stress tensor [32,38] of the polymer phase can be expressed as σ here, the Lamé coefficients, λ s and µ s and the elastic strain tensor of the solid phase, E are defined as [39] respectively; E 0 is Young's modulus and ν is the Poisson ratio.
The osmotic pressure inside the hydrogel comprises the mixing and ionic contributions, which can be expressed as where k B is Boltzmann's constant, V m is the equivalent volume occupied by one monomer, N x is the degree of polymerization, χ is the polymer-solvent interaction parameter, and c * k is the concentration in the fluid channel (outside hydrogels) of the ion species k. The polymer-solvent interaction parameter, χ(T, φ (p) ), is generally expressed as a function of temperature and polymer volume fraction [40][41][42]; i.e., where ∆h and ∆s denote the changes in enthalpy and entropy, and χ 2 is a parameter to express the polymer volume fraction dependence of the interaction parameter. The fixed charge concentration is given by [31] where c 0 f and c H + are the initial fixed charge and hydrogen ion concentrations, respectively, K a is the dissociation constant of the fixed charge group and H is the hydration state of the hydrogel, which is defined as the ratio of the volume of the fluid phase to the volume of the polymer phase inside the hydrogel, H = V f /V 0 . For axially restrained cylindrical hydrogels, hydration can be related to the strain of hydrogel as where E rr and E θθ are the radial and tangential strains, given by [43] E rr = ∂u r ∂r (17) respectively. Note that-compared to De et al. [31]-an additional term for hydration has been included: viz., tangential strain, since the deformation in the radial direction will trigger strain in the tangential direction [44]; hence, hydration should take into account total strain in radial and tangential direction. The dynamic viscosity of the fluid phase can be expressed as [35] µ (f ) = a 1 (T + a 2 ) a 3 (19) and the permeability of the polymer network is given by [45] where a i , κ (p) 0 and n are constants summarized in Table 1.
The effective diffusivity of ion inside the hydrogel is taken into account by Bruggeman equation [46]: where D k is the diffusive coefficient of ion species in water.
Here, pH and pKa are the negative logarithm of hydrogen ion concentration and dissociation constant given by respectively; c 1 is a constant presented in Table 1.

Boundary and Initial Conditions
The boundary conditions can be summarized as follows: • At the inlet of the channel, we prescribe • At the outlet of the channel, we prescribe • At the walls of the channel, we prescribe • In the centre of hydrogel, we prescribe u = 0 • At the hydrogels/fluid interface, the fluid velocity and fluid pressure are prescribed as Initial conditions invoked are Here, n is a unit vector normal to the given surface, | − and | + denote condition inside and outside the hydrogel, and c 2 is a constant presented in Table 1. The boundary condition in the centre of hydrogel is necessary in order to prevent translational movement of the hydrogel and corresponds to the way the hydrogel is attached to the flow channel; see e.g., [14,15].

Numerical Methodology
The mathematical model is solved with the commercial finite-element solver, Comsol Multiphysics 3.5a. Two geometries-hydrogels and channel-are solved simultaneously. Overall, the mathematical model for the hydrogels and flow inside the channel consist of eight dependent variables: u r , c 1 , c 2 , c 3 , u (f ) , v (f ) , w (f ) , and p (f ) . The geometries are resolved with around 1600-1800 elements to ensure mesh-independent solutions, amounting to around 2.6 × 10 4 -4.4 × 10 4 degrees of freedom; a finer mesh is chosen at the interface between a hydrogel and the surrounding fluid in the microchannel. The computations were carried out on a computer with a 2.66 GHz dual processor and 4 GB RAM and took around 10-30 min.

Calibration and Validation of the Hydrogel Model
Before we study the behavior of a hydrogel and its effect on the overall fluid flow in a microfluidic T-junction, we calibrate ∆h, ∆s, and χ 2 with the steady-state swelling curve for a diameter of 300 µm (training set) from the experiments by De et al. [31], as shown in Figure 2, and validate the deformation of the pH-sensitive hydrogel HEMA with experimental hydrogels with a diameter of 500 and 700 µm (test set). Overall, good agreement is achieved between the model prediction and the experiments. Clearly, the pHEMA hydrogel collapses at low pH and swells at high pH: This pH-induced swelling behavior can be attributed to the presence of acidic groups bound to the polymer chains, which become highly ionized at certain pH value [16,50]. The acidic group inside the hydrogel is only slightly ionized when the pH drops below the pKa of the hydrogel-in this case the pKa is 5. As a result, swelling of the hydrogel at pH changes below pKa is marginal, which is mirrored by the slight increase in the hydrogel radius for pH 4 in Figure 2. As the pH increases and approaches the pKa value, the acidic functional group becomes near to fully ionized by deprotonation, which results in an increase in the fixed charge density, c f . The increase in the fixed charge density, in turn, appears in the osmotic pressure (driving force) and causes a swelling, as depicted in Figure 2 for the pH range 4-6. When the ionization process reaches its saturation point, an increase in pH does not affect swelling behavior of the pHEMA hydrogel [51]. This can be observed in Figure 2 where the hydrogel stops swelling at leading order when pH > 7.  Turning our attention towards the deformation kinetics of the hydrogel, we first calibrate the permeability constant, κ 0 , for the shrinking of a 300 µm hydrogel when subjected to pH changes from 6 to 3 (training set), as illustrated in Figure 3a. We then validate the deformation kinetics with the corresponding swelling (see Figure 3b) with reasonably good agreement. Overall, we note that the shrinking is approximately ten times faster than the swelling: shrinking and swelling require around 1500 s and 18, 000 s respectively in order to reach the new steady state.
Aside from calibration and validation purposes, it is of interest to study equilibrium swelling behavior and deformation kinetics of a hydrogel since these are two important key factors in designing hydrogels for microfluidic flow control: From the equilibrium swelling behavior, we can identify and modify the properties of hydrogels that affect the swelling ratio and synthesize a hydrogel with the desired swelling ratio for flow control purposes; and from the deformation kinetic behavior, we can estimate the response time of the hydrogel when it is employed as microfluidic flow controller, after which we can design or synthesize a hydrogel with a sufficiently fast response for any given microfluidic flow control system.

Flow Behavior Inside a T-Junction with One or Several Hydrogels in One Branch
We proceed further by examining the sensing and actuating behavior of a 300 µm hydrogel in a T-junction when the solution pH is changed between 3 and 7, as illustrated in Figure 4. This configuration represents a simple microfluidic flow controller based on a stimuli-responsive hydrogel-also commonly referred to as resistance-based flow control [15]. Initially, at low pH, the hydrogel is in its shrunken state and thus allows fluid flow (Figure 4a) between itself and the walls of the microchannel. A step change in pH from 3 to 7 is then applied to the system, for which the hydrogel starts to swell towards the new equilibrium and block the channel, as depicted in Figure 4b-d. The mass flow rate of the fluid at the inlet and outlets of the channel as the hydrogel deforms is presented in Figure 5. When the hydrogels reach new equilibrium at t ∼ 200 min, a step change in pH from 7 to 3 is applied; thus, the hydrogel starts to shrink towards the initial condition. Here, several features are apparent: First, the mass flow rate at the inlet of the channel decrease as the hydrogel swells, similar to that at the left outlet; second, the response of the hydrogel is rather slow-it takes around 120 min for the hydrogel to fully close the channel-which can defeat the purpose of flow control. The reason for the first observation is simple: as the hydrogel swells, it obstructs the flow and creates high resistance for the fluid to flow because of its low permeability. The second observation suggests that we should modify the hydrogel microfluidic system to obtain higher response rates.
In light of the second observation, we demonstrate two strategies in improving a hydrogel's response rate: first, by replacing a single larger hydrogel with multiple smaller hydrogels, and second, by employing a hydrogel with higher permeability (macroporous hydrogel).
For the first strategy, we implement two hydrogels with the size of 150 µm and three hydrogels with the size of 100 µm to replace the 300 µm, as shown in Figure 6. In doing so, we find that the response time is approximately 3 times faster (for 150 µm hydrogels) and 6 times faster (for 100 µm hydrogels) compared to the corresponding case with a single 300 µm hydrogel, as depicted in Figure 5. The reason for this response time enhancement is the fact that by reducing the size of hydrogel, we shorten the diffusion path of the penetrating fluid, which, in turn, leads to a faster response by the hydrogel. It should be noted, however, that by reducing the size of the hydrogels, we may reduce the mechanical strength and stability, which are necessary for a microfluidic flow controller [16,27]; therefore, careful consideration has to be taken to ensure an optimum design.   The second strategy is achieved by implementing a 300 µm hydrogel with 10 and 100 times higher permeability than the base-case hydrogel, which can be realized by utilizing macroporous hydrogels; see, e.g., [53,54]. With this approach, we find that the response times are around 10 and 95 times faster compared to the flow control system with a lower permeability, as presented in Figure 7. The faster response can intuitively be explained by the fact that swelling and shrinking kinetics mainly depend on the permeability of the hydrogel: on one hand, a low permeability induces a high resistance to the penetrating fluid flowing into the hydrogels, which in turn result in slow deformation response; on the other hand, a hydrogel with a high permeability allows for easier fluid penetration [35]. The utilized macroporous hydrogels, however, should possess a sufficiently low permeability, because hydrogels with too high permeability might allow a non-negligible amount of fluid to flow through them, which would defeat the purpose of flow control.

Flow Behavior Inside a T-Junction with a Hydrogel in Each Branch
In this configuration, two hydrogels with 10 times higher permeability than the base-case hydrogel with opposite behavior are introduced in the microchannel: a positive pH-responsive hydrogel, which swells as the pH increases, and a negative pH-responsive hydrogel, which shrinks as the pH increases. The fluid flow is thus either directed to the left or the right channel depending on the pH, as illustrated in Figure 8. Recalling that the shrinking is approximately ten times faster than the swelling for the conditions and HEMA hydrogel considered in this study, we expect that the switching between the positive and negative response hydrogels will not be symmetric. This is indeed the case, as during the first few minutes, the left channel starts to open before the right channel is fully closed. This, in turn, leads to fluid flow through both branches and causes an increase in the mass flow rate at the inlet, as shown in Figure 9. Clearly, one has to be careful when designing a flow sorter with hydrogels since the latter may still allow fluid flow through an undesired channel during the first few minutes (depending on response rate).   As the pH-positive hydrogel swells further, the right branch starts to be blocked; therefore, we see a decrease in mass flow rate at the inlet. When the positive pH-sensitive hydrogel reaches a new equilibrium (t ∼25 min), the entire right branch has been blocked, forcing all the fluid through the left branch (Figure 8d), for which we observe that the mass flow rate at the inlet is equal to that of the left branch. A step change in pH from 7 to 3 is then applied for which the positive pH-sensitive hydrogel shrinks whereas the negative pH-sensitive hydrogels swells towards the new equilibrium conditions.
In this configuration, each hydrogel plays the role of a sensor, regulator and an actuator for the fluid flow commonly handled by three different components. This particular configuration could be implemented in chemical or biochemical applications, where precise pH control is required; for example, in sequence determination of protein and DNA analysis [21].

Conclusions
A mathematical model for hydrogels embedded in a microfluidic T-junction that takes into account conservation of mass, momentum and ions for laminar, incompressible flow and for the sensing/deformation of a pH-sensitive hydrogel has been derived and presented. The predicted swelling behavior of a pH-sensitive hydrogel was validated with steady-state and transient experiments and achieved good agreement. The model was then employed to study the deformation behavior of hydrogels at various pH values and their impact on the fluid flow inside the microchannel where they act as autonomous valves. Two configurations were considered: a T-junction with hydrogels in one branch and a T-junction with hydrogels in each branch. Overall, the model could provide an insight into the swelling/shrinking behavior of hydrogels, which act as autonomous microvalves at various pH values.
From the numerical investigation, it was found that the response rate of hydrogels subject to pH changes is slow, which could defeat the flow control purposes. As such, two strategies to improve the response rate of the hydrogels were proposed and demonstrated: First, by using smaller hydrogels and, second, by employing hydrogels with higher permeability. It was found that the response rate improved 9 times when the hydrogel's size was reduced to 100 µm from 300 µm, and it could be further improved (up to 95 times) when macroporous hydrogels with 100 times higher permeability were implemented. It should be noted, however, that smaller hydrogels tend to have weaker mechanical strength while macroporous hydrogel may allow fluid to flow through them, which would defeat the purpose of flow control. Therefore, careful consideration is required when designing and synthesizing hydrogels for microfluidic flow control applications.
Finally, we would like to highlight that the model is not limited to pH-sensitive hydrogels; it can be extended to other stimuli-responsive hydrogels such as thermo-, electric-, alcohol-, and glucose-sensitive hydrogels.   Figure 9. Response of the fluid flow during swelling and shrinking in a T-junction with two 300 µm hydrogels in each branch. The mass flow rate of the fluid for pH changes between 3 and 7 are at (−) the inlet, (−·) the left outlet, and (· · ·) the right outlet.