PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Ind Eng Chem Res. Author manuscript; available in PMC 2010 April 27.
Published in final edited form as:
Ind Eng Chem Res. 2008 December 17; 47(24): 10053–10063.
doi:  10.1021/ie070957b
PMCID: PMC2860336
NIHMSID: NIHMS86652

In Vivo Simulations of the Intravenous Dynamics of Submicron Particles of pH-Responsive Cationic Hydrogels in Diabetic Patients

Abstract

A mathematical model describing glucose-dependent pH swelling and insulin release is developed for pH-sensitive cationic hydrogels in which glucose oxidase and catalase have been immobilized and insulin imbibed. Glucose based swelling and insulin release are simulated for intravenously injected particles at various design conditions. The effects of particle size, the number of injected particles, insulin loading, enzyme loading, monomer functional group loading and pKa, and hydrogel crosslinking ratio on insulin release and glucose sensitivity are investigated in order to optimally design the device for use. Increased insulin infusion is shown to result from increasing the number of circulating gels, increasing the collapsed particle size, or by decreasing the crosslinking ratio of the system. Release duration is shown to be dependent only upon the particle size and the achievable diffusion coefficient of the system. Glucose sensitivity, as measured by gluconic acid production and by the system pH, are functions of glucose oxidase loading and the concentration and pKa of the monomer used in the hydrogel.

The necessarily submicron particle size results in very rapid device insulin depletion. When the device is designed without considering constraints, the resulting release profile resembles that of an on/off switching mechanism. Future work will focus on simulations of swelling and release when the device is implanted in an alternative administration site.

Keywords: Simulation, Hydrogels, Glucose, Insulin, Drug Delivery, Model

Introduction

The goal of implicit closed-loop glucose control in Type 1 diabetic patients is to be able to release a specific amount of insulin in response to the levels of glucose in the bloodstream without the need for equipment such as sensors or pumps. To this end, Peppas and collaborators [1, 2] proposed the use of cationic, pH-responsive hydrogels of poly(diethyl aminoethyl methacrylate-g-ethylene glycol monomethyl ether monomethacrylate) (poly(DEAEM-g-EGMMA)) as a means to provide such control.

A cationic hydrogel system in which glucose oxidase has been immobilized onto the polymer backbone via covalent bonding with acryloyl chloride would result in excess glucose in the environment being oxidized to form gluconic acid. The resulting acid would cause the local pH of the system to decrease, resulting in swelling of the hydrogel. This swelling results in an increase in the mesh size of the gel. Podual et al. has demonstrated both the ability to successfully immobilize glucose oxidase and catalase onto a hydrogel backbone [1] and the ability of such a gel to reversibly swell and deswell in response to a change in the glucose concentration [3].

When insulin is incorporated into the gel, the increased mesh size results in increased mass transfer of insulin from the gel. As glucose is removed from the environment, the amount of oxidation is decreased, resulting in a decrease in the gluconic acid concentration and an increase in pH. This results in reversible swelling and deswelling as a result of the concentration of glucose in the environment. Such a system allows insulin to be released entirely as a function of the glucose levels in the bloodstream, resulting in the desired closed-loop control.

Before such a device can be implemented, testing in patients would be required to verify efficacy in release. One way to reduce the amount of testing required is through the use of simulation. A mathematical representation of the hydrogel dynamics can be combined with a mathematical model of a Type 1 diabetic patient’s glucose and insulin dynamics to simulate the effects implementing the hydrogel system.

Such a hydrogel device model would have to describe the diffusion of glucose within the hydrogel, the oxidation of glucose to form gluconic acid, the elimination of gluconic acid via diffusion, the effect of gluconic acid on the local pH of the system, the effect of pH on the swelling of the gel, and finally the effect of swelling on the mass transfer coefficients of the various species involved. The approach taken here is to model the simplest dynamics first and then to add complexity in order to simulate the desired characteristics.

With a set of model equations developed, the next step in the modeling process is to estimate the model parameters in order to produce meaningful and accurate simulations. While actual parameters from experiments or literature would be ideal, engineering approximations are often necessary, along with sensitivity studies, in order assign model parameter values.

The first goal in modeling the gel dynamics is to verify that a lower limit mass transfer coefficient exists such that two primary features in insulin release result. First, in the absence of a disturbance change in glucose (i.e., a meal), the insulin should be released from the gel over a long period of time. Second, the release rate from a gel in the absence of a glucose disturbance should be at the physiologic levels of basal insulin release, either from a healthy pancreas or from a pump. Because it is not realistic to expect a constant infusion rate from a gel device, the more appropriate objective would that the gel should be able to demonstrate rapid oscillatory swelling and deswelling in order maintain glucose levels within a physiological basal range. This work aims to achieve such objectives.

First, a simple model is developed describing insulin and glucose dynamics within the gel in basal, high pH conditions. To mimic the action of a healthy pancreas, the authors have chosen to develop a model based on an intravenous system that releases insulin directly into the bloodstream. Second, an engineering analysis is performed using reported literature values and basic assumptions to determine reasonable parameter values for the simple model. Finally, the model parameters are varied by orders of magnitude to determine their effects on release with the minimal model [4], as they relate to the two desired primary features.

Methods

Model Development

To describe the dynamics of a feedback control system, we constructed the block diagram shown in Figure 1. When the glucose concentration in the blood is greater than normal levels, a concentration gradient exists between the glucose in the blood and the equilibrium levels within the polymer system. This results in the diffusion of glucose into the system. However, as glucose is diffusing within the gel, the presence of immobilized glucose oxidase results in the homogenous oxidation of glucose to form gluconic acid. This species exists in greater concentrations than that in the external environment. Tthe gluconic acid concentration in the blood is assumed to be zero, based on the idea that blood acts as a perfect sink for foreign hydrogen ions as the pH is rapidly regulated by the buffers in the blood. This results in passive diffusion of ionic species to the external environment.

Figure 1
Block diagram of the implicit closed-loop control process. Glucose diffuses into the hydrogel along a concentration gradient. Within the gel, glucose is oxidized to form gluconic acid, which causes the pH to decrease within the gel. The pH change results ...

The presence of gluconic acid results in a pH change at the local environment, resulting in a change in the size of the gel. This swelling also results in an increase in the diffusion coefficients of species through the hydrogel system, which in turn causes increased diffusion of insulin from the gel to the external environment. Therefore, in order to model the system, the diffusion of glucose, gluconic acid, and insulin must be modeled, along with the effect of gluconic acid on the local pH and the effect of pH on the size of the hydrogel.

Assuming no dynamics associated with metabolic glucose uptake or production, the dynamics of glucose in the external environment is described by Equation (1).

d(VBG)dt=A(NG)
(1)

Here, G is the blood glucose concentration, VB is the distribution volume of the bloodstream, and the mass transfer area, A, is given by Equation (2).

A=Ng4πR2
(2)

Here, Ng is the number of hydrogel spheres are in the bloodstream, and R represents the radius of a hydrogel particle. The term NG is the mass transfer flux from the bloodstream to the gel interior, as given by Equation (3).

NG=DGL(GGD)
(3)

In Equation (3), DG is the diffusion coefficient of glucose within the hydrogel network, L is the diffusion length across the gel, and GD is the glucose concentration within the hydrogel. For simplicity, spatial considerations have been eliminated and it is assumed that the internal glucose concentration beyond the characteristic diffusion length from the outside radius is uniform throughout the gel.

While a more accurate description would be one that included gradient terms, the purpose of the model is to investigate the feasibility of the hydrogel device as a control system, and a simplified model allows such feasibility to be determined. When the distribution volume of the bloodstream is assumed constant, Equation (1) through Equation (3) can be combined to give the following dynamic equation for glucose in the bloodstream:

dGdt=4πNgDGR2LVB(GGD)
(4)

As with Equation (1), the dynamics of glucose within the hydrogel can be described by a basic material balance, as shown in Equation (5).

d(VDGD)dt=A(NG)VDrG
(5)

In Equation (5), VD is the volume of the particle environment. The parameter rG is the rate of glucose oxidase-catalyzed glucose oxidation. This reaction is given [5] in Equation (a) and Equation (b).

Glucose+O2GlucoseOxidaseGluconicAcid+H2O2
(a)

2H2O2Catalase2H2O+O2
(b)

In the reaction, glucose is oxidized by glucose oxidase to form gluconic acid and hydrogen peroxide. As hydrogen peroxide is a toxic substance and must be eliminated from the body, catalase is also immobilized into the gel network. Catalase catalyzes the decomposition of hydrogen peroxide to form water and oxygen. Therefore, it is assumed that the oxygen content is enough to ensure that the glucose oxidation can occur freely, without oxygen limitations [1].

The rate of oxidation is described by Michaelis-Menten kinetics for enzyme catalyzed reactions and is given by Equation (6).

rG=VmaxGDKS+GD
(6)

Here, Vmax describes the maximum rate of reaction and is directly proportional to the amount of enzyme immobilized into the gel system. In addition, Ks describes the concentration of glucose when the reaction is proceeding at exactly half the maximum rate. The volume of the gel system, VD, is expressed by Equation (7).

VD=43πNgR3
(7)

Because the cationic hydrogel system will swell in response to pH, VD cannot be removed from the derivative.

Figure 2 shows the dynamic swelling response of two different hydrogel microspheres as described by Podual et al. [2]. When the particle size is on the order of 30 µm, the swelling and deswelling are showed to occur on the order of one minute or less. Since minute to minute control is the objective, and because the particles are an order of magnitude smaller than 30 µm, the assumption is made that the volume change will occur much quicker than the observed time scale of minutes. Therefore, a pseudo-steady state volume assumption can be made, allowing for the volume to be removed from the derivative. Combining Eq. (2) and Eq. (3) with Eq. (5), Eq. (6) and Eq. (7) results in the following simplified description for the dynamics of glucose in the hydrogel environment:

dGDdt=3DGRL(GGD)VmaxGDKS+GD
(8)
Figure 2
Dynamic swelling response of pH-sensitive hydrogel spherical particles to step changes in pH [2]. The circles represent a particle radius of 30 µm and squares represent a particle radius of 300 µm.

The dynamics of insulin within the gel is primarily driven by two processes. The first is the diffusion of solution insulin from the gel to the external environment. The second is the dissolution of any undissolved insulin in the hydrogel solution. If the assumption is made that the insulin concentration within the gel is below the saturation concentration of insulin (8 mg/mL), then the dynamics of insulin within the gel reduces to the following equation:

ddt(VDID)=A(NI)
(9)

The terms VD and A have previously been defined, ID is the insulin concentration in the gel, and NI is the insulin flux from the device, as described by Equation (10).

NI=DIL(IDI)
(10)

Using the simplifications of Equation (2) and Equation (7), and assuming again that the time scale of importance is greater than the time scale associated with the change in the particle radius, the dynamics of insulin within the gel can be simplified by combining Eq. (2), Eq. (7), and Eq. (9) to form the following expression:

dIDdt=3DIRL(IDI)
(11)

The insulin infusion rate from the gel is merely the absolute value of the dynamic change in insulin mass within the gel system, as described by Equation (9). Combining Equation (2), Equation (9), and Equation (10) results in the expression for the insulin infusion rate from the implicit closed-loop control system.

U=4πNgR2DIL(IDI)
(12)

Finally, the dynamics of insulin in the external environment can be determined in an analogous manner to the dynamics of glucose in the external environment. Neglecting the metabolic production and elimination of insulin within the body, the accumulation of insulin within the bloodstream depends entirely the rate of insulin infusion from the device.

ddt(VBI)=U
(13)

Assuming the circulation volume is constant, Equation (13) is simplified to the following dynamic expression for insulin in the bloodstream:

dIdt=UVB
(14)

It should be noted that the dynamics of insulin and glucose, as given by Equation (4) and Equation (14), assume no metabolic production or elimination of either species. In reality these effects clearly exist, and in order to use the dynamic equations to simulate the in vivo response of the system, a patient model must be chosen, and the form of the model chosen will determine the form of the glucose and insulin dynamic equations.

The ability of the gel to both sense glucose and change its transport properties in response to its presence is made possible by the production of gluconic acid. Gluconic acid dynamics was represented in a similar manner to the representation of glucose and insulin. The overall mass balance is described by Equation (15):

dGADdt=VmaxGDKs+GD+3DGARL(GADGA)
(15)

The gluconic acid production rate is equivalent to the device glucose oxidation rate. Gluconic acid is then assumed to passively diffuse down a concentration gradient into the external environment.

While the presence of gluconic acid is integral to the implicit control process, it is important to think in terms of pH for two reasons. First, experimental data shows the gel response as a function of pH. Second, and more importantly, the pH of the body cannot sustain significant pH changes, and therefore an analysis of the type of pH changes in the bloodstream must be performed in order to ensure safety during the device’s application.

There are two primary phenomena that determine the effect of the production of gluconic acid on the pH of the local environment of the gel system: the existence of biological buffers that keep the pH of the system tightly controlled, and Donnan equilibrium. The buffer system in equilibrium is described by the Henderson-Hasselbalch equation:

pH=pKa+log[BufferBase][BufferAcid]
(16)

When gluconic acid is produced in the gel, the gel buffer quenches the foreign acid to form more of the buffer acid, and the new pH is based on the new acid and base concentrations [5]:

pHnew=pKa+log[BufferBase]0GA[BufferAcid]0+GA
(17)

While the concentrations of the major buffer species in the blood are readily obtained [6, 7], these values are not the same inside the gel environment. Because there are negatively charged methacrylate groups immobilized on the polymer backbone of the gel, and because the system pH is greater than the pKa of the methacrylate backbone, there are different concentrations of ionic species in the gel compared to the bloodstream (Donnan effect [8]). When Donnan equilibrium is established, all charges are balanced at the cost of having different species concentrations across a semi-permeable membrane, which for this case is the gel device. The concentration differences for species across the membrane are described using the Donnan ratio, K [9]:

ci,gelci,blood=Kzi
(18)

For the ith ionic species, the ratio of gel concentration to the blood concentration is simply the Donnan ratio raised to the power of the valence of the ionic species. The Donnan ratio can be found by performing an overall charge balance to develop the following relationship:

(1φ)iziKzici,blood+σφ1+K110pHpKg=0
(19)

Equation (19) is based on the assumption of ideal Donnan equilibrium. As discussed by Siegel [10], the assumptions included ion activity equal to its concentration, that the fixed amine groups of the polymer backbone be assumed to have a single ionization constant, K (or pKg), that bulk electroneutrality does indeed hold, and that the species concentrations are given with respect to the fluid solvent within the gel, not the entire gel volume. The possibility of positive ions in plasma associating with the amine groups, resulting in the ideal Donnan equilibrium assumption being invalid, is also noted. However, provided there exists an anion in relatively high concentration, the assumptions are applicable. For the purpose of simulations, ideal Donnan equilibrium is assumed valid.

The variable pKg represents the pKa of the methacrylate functional group of the polymer backbone. The variable σ represents the density of functional groups in the polymer. This is a design parameter for this system. These two variables can be selected by choosing specific materials. The variable [var phi] represents the volume fraction of the polymer in the device system. Because the gel is assumed to be swollen as a result of absorbed water, the actual volume fraction of polymer will be significantly smaller than that of the collapsed particle. The volume fraction is represented by Equation (20):

φ=VpolymerVtotal=vQv=1Q
(20)

The volume fraction is simply the total polymer volume divided by the total volume. Given a collapsed polymer volume v, the total volume of the device is characterized by the volume swelling ratio Q, which is defined as the volume of swollen polymer complex to the initial volume of collapsed polymer:

Q=R3Rcollapsed3
(21)

The product Qv denotes the total device volume, and it is easily shown that the volume fraction is merely the inverse of the swelling ratio. Given value of swelling ratio, the polymer design characteristics, and the ionic concentrations of the different species in the blood and the gel, the Donnan ratio can be determined using a nonlinear system solver and used to determine the buffer acid and base concentrations within the gel, as well as the normal pH within the gel as a result. Given the gel pH and buffer concentrations, the pH as a function of gluconic acid concentration can then easily be determined.

Peppas and co-authors [1, 2] have demonstrated a dramatic increase in the swelling ratio for poly(DEAEM-g-EGMMA) when subjected to a pH below 7.1. These results are given in Figure 3. Using the pseudo-steady state assumption with respect to volume swelling, it was assumed that the swelling ratio will rapidly approach its equilibrium value as shown in Figure 3. For the purposes of simulation, only empirical models were needed to describe the effects of pH on swelling. Based on the form of the curve of the results shown in Figure 3, a hyperbolic tangent relationship was used, as indicated below:

Q=A+Btanh(C·pH+D)
(22)
Figure 3
Equilibrium swelling response of pH-sensitive cationic hydrogels as a function of pH and crosslinking ratio [1]. Circles represent a crosslinking ratio of 0.005, squares represent 0.01, diamonds represent 0.02, and triangles represent a ratio of 0.04. ...

It can be seen in Figure 3 that a hyperbolic tangent relationship does not perfectly match the swelling profile, as the gel continues to swell as the pH drops. However, because the pH of the system will not deviate from the range of 7.0 in the gel to 7.4 of the blood, it was concluded that such a relationship is adequate for the interval selected.

Figure 3 also shows the swelling effects of pH for gels of different crosslinking ratios. As the crosslinking ratio is increased, the device becomes much more stiff. To quantify this effect, empirical relationships using the data of Figure 3 were developed. Based on the hyperbolic relationships, at maximum and minimum swelling, the hyperbolic tangent function can be approximated as unity. This allows A and B to be determined based on the asymptotic swelling values.

A+B=Qmax
(23)

AB=Qmin
(24)

From an inspection of Figure 3, Qmin does not change significantly as a function of crosslinking ratio, and it was approximated as 3 at crosslinking ratios lower than 0.02 and 2 at crosslinking ratios of 0.02 or greater.

Clearly Qmax is strongly dependent on crosslinking, and the maximum swelling values were plotted as a function of crosslinking ratio, as shown in Figure 4. A decaying exponential was chosen to describe the data, simply because it gave the best fit. The relationship is described by Equation (25).

Qmax=40exp[49.8(X0.005)]
(25)
Figure 4
Maximum swelling ratio of pH-sensitive cationic hydrogel as a function of the crosslinking ratio. The line represents an exponential function used to fit the data.

This fit was based on the data observation that a crosslinking ratio, X, of 0.005 (moles crosslinking agent/moles monomer) has a maximum swelling ratio of approximately 40. As shown in Figure 5, parameter C exhibited no crosslinking dependence, and so an average value of −10 was chosen. Finally, because the swelling transition is supposed to occur at the gel pKg, the hyperbolic tangent functions inflection point was assumed to occur when the pH is equal to the pKg. Taking the second derivative of the function shows that the hyperbolic tangent function goes to zero at this point. To show this mathematically, the parameter D is described by Equation (26).

D=C·pKg
(26)
Figure 5
Change of hyperbolic tangent term C of Eq. (22) with respect to crosslinking ratio, X.

With parameters A, B, C, and D expressed as a function of crosslinking ratio, the swelling ratio for any pH can now be calculated as long as the crosslinking ratio is known. Given the swelling and crosslinking ratio, the remaining transport parameters can be determined as well. Diffusion of a species within the gel is based on the size of species relative to the size of the mesh of the network. The mesh is described by Equation (27) [11].

ξCn1/2Q1/3N1/2lc
(27)

Equation (27) is based on the assumption of isotropic stretching of the polymer chains during the swelling process. The term Q1/3 represents the relative increase in size of a chain. The parameter Cn is known as the characteristic ratio and is taken to be 14.4 for poly(DEAEM-g-EGMMA) [1]. The term lc is simply the length of the bond between two carbon atoms (1.54 A). The term N is the number of repeating units between two crosslinks. Combined, N and lc give the average length of the unstretched chain. Each cross-linking agent has four functional groups, enabling to link four polymer chains together. However, each chain is shared between two crosslinking agent functional group sites. Thus, N can be approximated as shown:

N12X
(28)

Given the mesh size, the swelling ratio, and the size of the particle being transported, the diffusion coefficient for a species i through the gel is given by Equation (29) [11].

DiDi,0(1rH,iξ)exp(YQ1)
(29)

Y is a constant assumed to be one for this and most other polymeric systems [1], and rH,i is the hydrodynamic radius of the particular species. Equation (29) is a scaling law based on statistical thermodynamics and the probability that both enough free volume exists for the solute to move through the gel and that the gel network conformation is such that the solute can move through it as well.

The final process associated with model development is with respect to the degradation of the polymer. As an approximation of degradation, the works of Batycky et al. [12] and Martens et al. [13] have been studied and adapted for crosslinked gels. Batycky et al. showed that microparticles of 50:50 poly(D,L-lactic-co-glycolic acid) (PLGA) degrade by the cleaving of bonds on the polymer chains, resulting an increase in the porosity of the spheres. This bond cleavage was modeled as a zero order process because all bonds are equally likely to be cleaved. Martens and co-authors investigated the bulk degradation of hydrogels of poly(lactic acid-b-ethylene glycol-b-lactic acid) (poly(LA-b-EG-b-LA)) and assumed first order kinetics with respect to the degradable units.

To develop pH-responsive hydrogels that are degradable, a degradable crosslinking agent must be used. As the system degrades, instead of the average pore size increasing, the mesh size will increase. With respect to bond cleavage, the only degradation that occurs will be respect to the bonds between a chain and the crosslinking agent. Based on this analysis, the degradation rate was assumed to be first order in the crosslinking ratio only.

dX(t)dt=kdX(t)
(30)

As the ratio decreases, insulin will be able to diffuse faster from the gel, until a certain threshold value is reached in which the insulin will be rapidly dumped from the device. The goal of design will then be to ensure that insulin is essentially depleted from the device before such degradation occurs.

Finally, given the model of hydrogel dynamics, implicit closed-loop control can be investigated by incorporating the model into a patient model. For this purpose, the minimal model [4, 14] was used.

dG(t)dt=P1(G(t)Gb)G(t)X(t)+D(t)
(31)

dX(t)dt=P2(X(t)Xb)+P3(I(t)Ib)
(32)

dI(t)dt=nI(t)+U(t)V
(33)

Parameter Determination

To perform the simulations, numerical values are needed for each parameter. Many of the parameters could not be determined from experiment, requiring other methods. Model parameters were determined in three ways. First, available values from literature were used when applicable. Second, parameters from the simulation results [5, 15] were used as first approximations for some of the values. Finally, intuitive reasoning was used to determine approximate values of the remaining parameters.

Furler et al. [14] provided the basal concentrations of glucose and insulin, as well as the parameter values of the patient model. Other research [5, 15] was used to determine both the hydrodynamic radius for the different species and their associated solvent diffusion coefficients, Di,0. The degradation rate, kd, was approximated as the same order of magnitude as the random scission rate given in reference [12]. The parameter values for Vmax, Ks, pKg, and σ were initially approximated using the parameters of reference [5]. Finally, the concentrations of buffer species could be determined [6, 7].

To determine the maximum loading of insulin within the gel, the saturation concentration of insulin in water was used [15]. It was then assumed that insulin is always dissolved within the gel. The average particle size was given an upper limit of 5 µm. In order to be able to circulate in the bloodstream, the particles must be able to flow through the capillaries. Given the average capillary internal diameter of 5–9 µm [6], the gel would likely have to be much smaller than this value. The nominal values of all parameters are summarized in Table 1.

Table 1
Nominal Parameter Values Used For Hydrogel Device Simulations

The final parameter, Ng, was determined based on the particle size and some assumptions about the administration of the device. It was assumed that the gels would be injected into the bloodstream in a solution that does not cause them to swell initially. A patient using insulin injections will likely inject between 3/10 to 1 mL of U-100 insulin daily for basal insulin [16]. U-100 means that the insulin solution contains 100 units per mL of solution. A device designed to provide insulin for multiple days would have inject a higher volume of solution every few days. As an initial assumption, it is assumed that 10 mL of solution will be injected into the bloodstream periodically. To be comparable to intravenous injections of other species, a basis of glucose [17], which has previously been injected in solutions at up to 50 wt %, was chosen. To be consistent with glucose, it was assumed that the particles should only represent up to 50 wt% of the injection solution. Assuming the particles have a density of approximately unity, the entire volume of particles can only equal 5 cm3. Therefore, Ng is determined by noting that the total number of particles must have a volume no greater than 5 cm3. Table 2 shows the maximum number of assumed particles for different particle sizes.

Table 2
Particle Size Effects on Maximum Number of Circulating Particles

It should be noted that, assuming rapid diffusion and rapid oxidation, the gluconic acid concentration within the gel would reach a maximum value equal to the glucose value in the blood stream. In the worst case, a scenario would occur in which all the gluconic acid from the gels is released into the bloodstream. Given that the total particle volume is only 5 cm3, compared to the body’s volume of more than 3 L of plasma [8], the acid concentration would likely be considerably diluted. Finally, considering the plasma carbonate buffer concentration of 25 mmol/L, it can be assumed that gluconic acid would be quenched nearly instantaneously with nearly no impact on the buffer system. Therefore, the gluconic acid concentration in the bloodstream is considered to be zero at all times.

Finally, the model neglects any immune response effects on the delivery system. The presence of ethylene glycol grafts has been assumed to provide stealth characteristics, although the duration of the system in blood has not been determined. The small size of the particles also means that they could be eliminated via opsonization [18], in which the particles are enveloped and ultimately eliminated by the body. While these concepts are extremely important with respect to the system design, it is assumed here, for the sake of investigating feasibility with respect to insulin release and glucose control effectiveness, that immune response effects are not present.

Simulation Objectives

The objectives of model simulation were to investigate the effects of each parameter on the system and to design a gel that is able to effectively provide adequate glucose control. To achieve satisfactory control, three criteria were used. First, insulin release from the gel should be sustainable for at least two to three days. A regimen requirement daily injections is no improvement over the currently used regimen. Second, the device should be able to provide a basal level of insulin throughout the day. As a specified basal value cannot be expected from a polymeric device, this criterion establishes that glucose should stay within a normal range of 70 to 100 mg/dL (3.9 to 5.6 mmol/L) during the day. Finally, the response of the device to a glucose infusion should result in a rapid increase in insulin diffusion, in order to provide adequate control of the disturbance. Using the simulation results and the control criteria, the parameters were adjusted in order to select the optimal gel design that results in effective control.

Computational Methods

All simulations were performed in MATLAB. To solve the differential equations, the MATLAB function ode15s was used, which is based on a modified Rosenbrock function to solve stiff systems [19]. To solve the system for the Donnan ratio, the function fsolve in the MATLAB Optimization Toolbox was used.

An investigation of the system equations shows that swelling is a function of pH. However, the pH is a function of the buffer concentration in the gel, and the buffer concentration is a function of the Donnan ratio. Because the Donnan ratio is a function of the swelling ratio, the system had to be solved iteratively in order to simultaneously determine the gel pH and the swelling ratio.

Finally, the system was initialized by providing the design characteristics and the basal values of the patient model. The initial values of the device insulin and crosslinking ratio were provided. The initial concentrations of device glucose and gluconic acid were determined using iterations. Because transport was pH dependent, the initial pH, gel swelling, diffusion coefficients, Donnan ratio, and state values had to be determined simultaneously.

Results

Parameter Effects on Device Insulin Release

Before the gel’s ability to control glucose levels is investigated, the system must be designed such that insulin is able to be infused at realistic levels for a period of multiple days. Thus, the transport of insulin was investigated initially to ensure that sufficient release can be achieved. The primary design characteristics important to release are the duration of the release from the device and the magnitude of the insulin infusion rates. All investigations were performed at basal conditions.

Figure 6 shows the effects of insulin loading on the device behavior. The upper limit on insulin loading was established to be the saturation concentration of insulin in water. Figure 6 shows that insulin loading had no effect on the duration of release. Insulin loading did affect the flow rates that were achievable. Mathematically, this is easily seen in Equation (12), where the concentration driving force is directly dependent on the insulin concentration within the gel. Increasing or decreasing the insulin loading causes a directly proportional increase or decrease in the release rates.

Figure 6
Fraction of insulin released and insulin infusion rate as functions of time and the initial gel insulin concentration.

Figure 7 shows the effects of the number of particles assumed to be circulating at any time. Because the total number of particles possible is size dependent, these results were all determined for a gel with an average collapsed radius of 1.5 µm. As with insulin loading, the number of gels circulating played no role in the duration of release. This can be explained simply by investigating Equation (11). Assuming uniform insulin concentrations, each gel will have the same insulin profile, dependent only upon the mass transfer coefficient and the concentration gradient of insulin. Also like insulin loading, the number of gels has a directly proportional effect with respect to changing the infusion rate from the polymeric device. This can also be seen in Equation (12), in which insulin infusion is directly proportional to the number of circulating particles.

Figure 7
Fraction of insulin released and insulin infusion rate as functions of time and the number of circulating gels. The particle had a collapsed state radius of 1.5 µm.

Figure 8 shows the effects of the collapsed particle radius on the gel release characteristics. For each simulation, the number of gel particles in circulation was set at 3.54×107. The maximum Ng for a particle of radius 1.5 µm is 3.54×1010, and the number was decreased by 3 orders of magnitude in order to ensure the number was feasible for bigger particles, in accordance with Table 2. The figure shows a strong dependence on particle size. As the particle size was increases, the duration of insulin release increased. This was a result of the increased diffusion path from gel to the circulation. Mathematically, it is seen in Equation (11) that the gel insulin concentration change is inversely proportional to the particle size. The infusion rate, however, will depend on the flux and the gel surface area, resulting in a directly proportional relationship between insulin infusion and the collapsed size of the particle.

Figure 8
Fraction of insulin released and insulin infusion rate as functions of time and collapsed particle size. For all sizes, there were 3.54 × 107 particles in circulation.

Figure 9 shows the effects of the diffusion coefficient on the gel characteristics. The duration of release was a strong function of the order of magnitude of the diffusion coefficient. An order of magnitude decrease in the diffusion coefficient resulted in an increase in the duration of release from the microparticle. As shown in the figure, decreasing the diffusivity also resulted in a decrease in the insulin infusion rate, although the device was able to infuse this level of insulin for a longer period of time. Mathematically, this can be analyzed by an investigation of the model equations. As insulin flux is proportional to the diffusion coefficient, a higher coefficient will result in a higher insulin release rate. The higher infusion rate, when the particle size, the number of gels, and the insulin loading are maintained constant, means the gels are losing their insulin at a faster rate, thus leading to a shorter duration of release.

Figure 9
Fraction of insulin released and insulin infusion rate as a function of time and the diffusion coefficient of insulin, Di,0.

While the diffusion coefficient of glucose could also have been investigated, it was assumed here that it will likely always be approximately an order of magnitude larger than the insulin coefficient, as a result of the much smaller hydrodynamic radius and its known diffusion coefficient in water relative to that of insulin. Thus, a change in the glucose coefficient was implied when the change was made in the insulin coefficient.

Figure 10 shows the effects of the gel crosslinking ratio on the gel insulin release characteristics. As the figure shows, decreasing the initial crosslinking ratio resulted in an increase in the insulin infusion rate and a decrease in the duration of the gel release, as more was released at a time. Mathematically, this can be explained by observation of Equation (27) and Equation (28). Physically, as the crosslinking ratio decreases, the average gel mesh size will be larger, as there will be more repeat units between crosslinks. In addition, as previously mentioned, the lower crosslinking ratio will also lead to higher swelling ratios for the same pH environment. Because the swelling ratio also contributes to the mesh size, the crosslinking ratio affects the mesh size in both fashions. An investigation of Equation (29) shows that both a larger mesh size and a higher swelling ratio lead to an increase in the diffusion coefficient.

Figure 10
Fraction of insulin released and insulin infusion rate as functions of time and the initial crosslinking ratio of the hydrogel network. The rapid fall is a result of the high infusion rates causing the particle to collapse, thereby decreasing the infusion ...

An interesting phenomenon was the sharp decrease in the insulin infusion rate in Figure 10. This occurred because of the internal feedback relationship between the pH and the swelling ratio of the polymer system. The pH determined the swelling ratio of the species. However, because the swelling ratio also contributes to the Donnan ratio, the swelling ratio affects the internal pH of the system. Thus the pH change causes a sharp drop in the swelling ratio, which contributes to further increasing the pH. Because the infusion rate is directly related to the particle size, the rapid decrease in the swelling ratio results in a rapid decrease in the infusion rate. It can also be seen in the case of X = 0.004 that as the insulin infusion rate decreases below the steady state requirement for the system, glucose levels increase, resulting in a decreased pH, eventually resulting in rapid swelling, and a rapid increase in the infusion rate. This oscillatory behavior allows the system to adjust its infusion rate in response to the environmental conditions.

Parameter Effects on Glucose-Induced Swelling

The parameters investigated were those primarily responsible for the insulin release magnitude and duration. In addition, to provide minute-to-minute control of glucose, the gel must be able to rapidly swell and de-swell in response to changing glucose conditions. While the diffusion of glucose is important, this transport will be related to the transport of insulin, and the gel will have to be designed based on the glucose diffusion rate that allows insulin to diffuse at a reasonable level. Of more importance are the effects of glucose on swelling and release. Given that release will be based on the previously mentioned parameters for a given swelling ratio, glucose sensitivity is primarily based on swelling. While realistically, glucose sensitivity would also depend on thermodynamic interactions with other components of the system, the model neglects these interactions. Thus, the important design characteristics are the rate of glucose oxidation to form gluconic acid and the resulting pH from the formation.

Figure 11 shows the effect of glucose oxidase immobilization with respect to gluconic acid formation and pH. As the figure shows, gels with enzyme loadings greater than 0.01 µM did not display different glucose sensitive characteristics by increasing the enzyme loading. At this loading upper limit, all glucose was oxidized very rapidly. As the figure shows, high loading resulted in a high gluconic acid concentration, whereas low loading resulted in a low acid concentration. This resulted in a low pH at basal conditions for high loading, and a high pH at basal conditions for low loading. For low loading, the high pH resulted in low insulin release, which caused glucose levels to rise. However, the low loading resulted in only a small amount of gluconic acid produced for the glucose rise, and the pH decrease was very slow. For high loading, the insulin infusion rate was initially high, causing glucose to fall. The drop in glucose caused a proportional drop in gluconic acid, raising the pH rapidly. This rise causes insulin release to decrease, resulting in a rise in glucose. This process cycles between high and low release rates, with the goal being to keep glucose oscillating in a physiologically acceptable range.

Figure 11
Gluconic acid concentration and Gel pH as functions of time and enzyme loading. The initial value is the gluconic acid concentration at basal conditions. For low loading, low gluconic acid levels result in low infusion rates, causing a rise in glucose ...

Figure 12 shows the effects of varying the pKg of the monomer used in the hydrogel. As can be seen, the pH in the device was below the transition pH of all three devices at basal conditions. This resulted in a large insulin infusion rate, causing glucose, and ultimately gluconic acid, to decrease at a rapid rate. This decrease caused the pH to increase. As the transition pH was reached, the system collapsed, and gluconic acid ultimately began to increase. As the figure shows, however, when the pKg of the monomer increased, it became more difficult to ever approach the transition pH. Given the conditions, the monomer of pKg 7.5 was not able to approach the transition pH value before gluconic acid was completely depleted. Therefore, no transition would ever occur, and this system would not be able to collapse.

Figure 12
Gluconic acid concentration and gel pH as functions of time and the pK of the monomer used in the hydrogel device.

Figure 13 shows the effects of the functional group concentration on the glucose response. Given the same initial gluconic acid concentrations, an increase in the functional group loading resulted in an increase in the basal state pH of the system. This allowed the state of the gel to be controlled such that the device was in the collapsed state at basal conditions. Mathematically, this is because the functional group loading directly affects the Donnan ratio, which is used to determine the buffer species concentrations within the gel.

Figure 13
Gluconic acid and gel pH as functions of time and the functional group concentration in the hydrogel device.

Optimization of the Hydrogel Device for Insulin Release

The optimal design of the device occurs in three phases. The first phase is to ensure that release can happen for a long enough duration. The second phase is to then make sure that the steady release is occurring at a high enough delivery rate to be physiologically beneficial. Finally, the system must be optimized in order to ensure glucose response occurs rapidly and effectively enough to maintain normal levels during basal conditions.

The duration of insulin release was shown to depend only upon the diffusion coefficient and the particle size. For an intravenous device, the particle size is limited by the size of the capillaries. Given an average capillary diameter of between 5 and 9 µm, the gels have to be able to flow through the capillaries without getting stuck. Realistically, the particles should then be designed on the submicron level. Using 5 µm as the upper limit for the swollen device diameter results in the maximum collapsed particle radius of 1.35 µm. However, such a device would only be able to pass through capillaries one gel at a time, and as such the gels would likely not be circulating at a uniform number concentration, as assumed for a pharmacokinetic model.

Based on the development of Equation (29), Di,0 is actually the diffusion coefficient in the pure solvent. As shown in Table 2, this value is on the order of 10−7 cm2/s. As shown in Figure 9, this is four orders of magnitude greater than a coefficient value that causes a larger particle to be depleted within two hours of circulation. Insulin release will not be at the desired duration for diffusion coefficients on the order of 10−15 cm2/s, a value that is not feasible for a mesh particle system. On the other side, assuming the diffusion coefficient in pure solvent was used, the particles would have to be on the order of one mm in diameter. Such a system would not be able to circulate. In addition, the pseudo-steady state assumption would no longer be valid, and the time constants and delays associated with viscoelastic swelling and relaxation would have to be considered. Thus, based on the transport requirements, such a device is not feasible for use as an intravenous delivery system.

While the previous paragraph essentially renders the hydrogel system infeasible as an intravenous insulin delivery device, optimization continued with the assumption that a gel could be developed with the necessary diffusion coefficient-particle size combination that results in release for a duration of three days or longer. Given the release duration, the necessary insulin infusion rate could be achieved simply by increasing or decreasing the insulin loading or number of particles accordingly. For the case in which large particles and fast diffusion is assumed, insulin loading and the number of particles will have to be reduced by several orders of magnitude to achieve desired levels, as infusion rates will be very high. In contrast, the slow diffusion and small particle assumption will require maximum insulin loading and the maximum number of gels that can feasibly be in circulation.

To optimize the gel for glucose response, two primary objectives must be achieved. First, the gel must be collapsed at basal levels of glucose and lower, and it must expand upon achieving high glucose levels. Second, this response must be fast. To achieve the second objective requires a high enough enzyme loading to result in a gluconic acid concentration within the gel that is able to cause pH changes in the presence of buffers. The first objective can be achieved by careful selection of the type and amount of monomer used in the device. For this optimization, the desired swelling ratio of 3 in the basal state was specified.

Based on this value, different combinations of pKg and functional group concentration were used in the determination of the Donnan ratio and ultimately the pH. This pH was then used to determine the swelling ratio of the system. If this value did not match the original value, the process was used again until the swelling ratios converged. As was often the case, this iterative process diverged to a large swelling ratio, resulting in low pH values. Figure 14 shows a plot of the specified swelling ratio vs. the actual ratio for a functional group concentration of 1 M. The large slopes usually resulted in divergence to the flatter regions of the graph.

Figure 14
Results of device optimization for optimal pKg and functional group concentration. For each pKg value, a Q is specified and used with a functional group concentration of 1 M to determine the Donnan ratio. The Donnan ratio is then used to calculate pH ...

Figure 15 shows the same iterative routine used with a monomer loading of 4 M. At the low swelling ratio region, the curves became flatter, indicating that the higher monomer loading allows the collapsed state to be achieved at the basal condition. Based on the figure, the first glucose response objective was achieved by using a monomer with pKg of 7.0–7.1, at a loading of 4 M or greater. Based on this optimization sequence, and the assumption that size and diffusion limitations are neglected, the gel insulin infusion and glucose response are shown in Figure 16 for basal conditions.

Figure 15
Results of device optimization for optimal pKg and functional group concentration. For each pKg value, a Q is specified and used with a functional group concentration of 4 M to determine the Donnan ratio. The Donnan ratio is then used to calculate pH ...
Figure 16
Minimal model implicit closed-loop controlled glucose response at basal conditions, over a period of 2.8 days. NG = 5×1010, DI,0 = 2.82×10–15 cm2/s, Rcollapsed = 937 nm, σ = 4 M, enzyme loading = 10 µM, and pKg ...

The change in glucose resulted in an order of magnitude change in the insulin infusion rate. In this way, the control system was most similar to a simple switching mechanism as opposed to a real controller. Thus, the implicit closed-loop control system was quite inferior to a controller developed using process control principles. Furthermore, the oscillations that occurred resulted in alternating conditions of hyperglycemia and hypoglycemia, although this behavior may have been a consequence of using the minimal model.

Conclusions

A model has been developed describing the processes associated with implicit closed-loop control, including glucose diffusion into the gel, the oxidation of glucose to form gluconic acid, the effect of gluconic acid on the gel pH, the gel pH effect on gel swelling, and swelling effects on transport. Each model parameter was varied by several orders of magnitude in order to show parametric effects on release duration, insulin delivery rates, and glucose responsiveness. Given the fundamental parametric understanding, the gel was then optimized for glucose control in Type I diabetic patients. Based on physiological considerations, the requirement of small particles and the inherently fast diffusion times associated with a mesh network result in the conclusion that the use of pH-responsive hydrogels as an intravenous implicit closed-loop glucose control system is likely infeasible, as delivery cannot be sustained for a long enough duration to offer any advantage to current treatment methods.

Acknowledgements

This work was supported in part by the National Institutes of Health (grant EB 00146 to NAP) and a National Science Foundation Fellowship to TGF.

References

1. Podual K, Doyle FJ, III, Peppas NA. Preparation and Dynamic Response of Cationic Copolymer Hydrogels Containing Glucose Oxidase. Polymer. 2000;41:3975.
2. Podual K, Doyle FJ, III, Peppas NA. Dynamic Behavior of Glucose Oxidase-Containing Microparticles of Poly(Ethylene Glycol)-Grafted Cationic Hydrogels In an Environment of Changing pH. Biomaterials. 2000;21:1439. [PubMed]
3. Podual K, Doyle FJ, III, Peppas NA. “Glucose-Sensitivity of Glucose Oxidase-Containing Cationic Polymer Hydrogels Having Poly(Ethylene Glycol) Grafts,” J. Control. Release. 2000;67:9–17. [PubMed]
4. Bergman RN, Ider YZ, Bowden CR, Cobelli C. Quantitative Estimation of Insulin Sensitivity. Am. J. Physiol. 1979;236:E667. [PubMed]
5. Abdekhodaie MJ, Wu XY. Modeling of a Cationic Glucose-Sensitive Membrane With Consideration of Oxygen Limitation. J. Membrane Sci. 2005;254:119.
6. Guyton A, Hall J. Textbook of Medical Physiology. 11th Ed. Philadelphia, PA: Elsevier Saunders; 2006.
7. Kaplan LA, Pesce JA, Kazmierczak SC, editors. Clinical Chemistry: Theory, Analysis, Correlation. 3rd Ed. St. Louis, MO: Mosby; 1996.
8. Cooney DO. Biomedical Engineering Principles: An Introduction To Fluid, Heat, and Mass Transport Processes. NY: Dekker; 1976.
9. Ricka J, Tanaka T. Swelling of Ionic Gels: Quantitative Performance of the Donnan Theory. Macromolecules. 1984;17:2916.
10. Siegel RA. Hydrophobic Weak Polyelectrolyte Gels: Studies of Swelling Equilibria and Kinetics. In: Dusek K, editor. Advances in Polymer Science 109, Responsive Gels: Volume Transitions I. NY: Springer-Verlag; 1993. pp. 233–267.
11. Lustig SR, Peppas NA. Solute Diffusion in Swollen Membranes. IX. Scaling Laws For Solute Diffusion In Gels. J. Appl. Polym. Sci. 1988;36:735.
12. Batycky RP, Hanes J, Langer R, Edwards DA. A Theoretical Model of Erosion and Macromolecular Drug Release From Biodegrading Microspheres. J. Pharm. Sci. 1997;86:1464. [PubMed]
13. Martens P, Metters AT, Anseth KS, Bowman CN. A Generalized Bulk-Degradation Model For Hydrogel Networks Formed From Multivinyl Crosslinking Molecules. J. Phys. Chem. B. 2001;105:5131.
14. Furler SM, Kraegen EW, Smallwood RH, Chisolm DJ. Blood Glucose Control by Intermittent Loop Closure in the Basal Mode: Computer Simulation Studies with a Diabetic Model. Diabetes Care. 1985;8:553. [PubMed]
15. Podual K, Doyle F, III, Peppas NA. Modeling of Water Transport In and Release From Glucose-Sensitive Swelling-Controlled Release Systems Based On Poly(Diethylaminoethyl Methacrylate-g-Ethylene Glycol) Ind. Eng. Chem. Res. 2004;43:7500.
16. Becton-Dickinson Syringe Capacity and Dose Scale. [accessed 6/18/07]. http://www.bddiabetes.com/us/main.aspx?cat=1&id=252.
17. Lozner EL, Winkler EW, Taylor FHL, Peters JP. The Intravenous Glucose Tolerance Test. J. Clin. Invest. 1941;20:507. [PMC free article] [PubMed]
18. Owens DE, III, Peppas NA. Opsonization, Biodistribution, and Pharmacokinetics of Polymeric Nanoparticles. Intern. J. Pharm. 2006;307:93. [PubMed]
19. Shampine LF, Reichelt MW. The MATLAB ODE Suite. SIAM J. Sci. Comput. 1997;18:1.