Swellable Hydrogel-based Systems for Controlled Drug Delivery

Abstract

The controlled delivery of drugs can be effectively obtained using systems based on hydrogels. Tablets, to be orally administered, represent the simplest and the most traditional dosage systems based on hydrogel. Their formulation and preparation require to mix and to compress, in proper ratios, various excipients, including a swellable polymer and a drug. Carriers for controlled release systems are usually cross-linked polymers able to form hydrogels that show peculiar release mechanisms, where both diffusion and tablet swelling play important roles.

1. Introduction

Hydrogels are hydrophilic polymer networks, able to absorb large amounts of water, increasing their volume (i.e., they are able “to swell,” giving rise to a phenomenon known as swelling) [1-3]. Networks can be composed of homopolymers or copolymers, and their network structure and physical integrity are due to the presence of cross-links, of chemical (tie-points, junctions) or physical (entanglements, crystallites) nature. Based on the stability of these cross-links, the hydrogel can withstand exposure to water; otherwise, it can degrade and dissolve in water, after a given exposure time.

The first known mention to potential use of hydrogels in medical science is probably due to Wichterle and Lim in 1960 [4], which proposed the use of methacrylate polymers for biomedical applications (i.e., the filling of eyes after enucleation). Since then, countless applications of hydrogels in pharmacy and medicine have been proposed, tested, and in several cases implemented as in common uses. Of course, during these (more than) 50 years, several very valuable reviews and reference books have been published on the work done too. A comprehensive list of these reviews/book today could be really huge, and it is beyond the scope of this work. Just a few examples are given in the following. With reference to the medical applications of hydrogels, the reader can refer to the handbook by Peppas [5], as well as to the reviews by Hoffman [1], by Peppas et al. [3], and by Cabral and Moratti [6]. With reference to the applications of hydrogels in drug delivery, the reader can refer to the reviews by Kamath and Park [7], by Peppas [8], by Peppas et al. [2], and by Hoare and Kohane [9]. It is worth to note that the scientific interest on the hydrogels and their applications is really large, witnessed not only by the large amount of papers published on the topic but also by the fact that several papers related to hydrogels are “most cited” in bibliographic databases. In 2012, the journal Advanced in Drug Delivery Reviews has devoted a special issue to the most cited papers (over than 1 thousand citations each) published on the journal in its 25-year lifetime, and four of the 32 papers were specifically devoted to hydrogels [1, 10-12].

1.1. Swelling hydrogels: Nature and classification

Hydrogels are interesting candidates for pharmaceutical and biomedical applications because of their characteristics such as:

1. Biocompatibility and bioabsorbable nature

2. Mechanical properties

3. Degradation properties

4. Compliance with sterilization protocols [7]

In particular, the balance between mechanical and degradation properties plays a relevant role in the design of a pharmaceuticals as well as a scaffold made of hydrogel to allow the desired release profile from the pharmaceutical in the time interval needed by the drug to perform the desired effect on the body or to ensure the structural function of the scaffold over the necessary time interval for the body to reconstruct the damaged part.

Several hydrophilic polymers can be used to produce hydrogels [1]:

1. Natural polymers, for example, alginate and carrageenan (anionic polymers); chitosan (cationic polymer); dextran, agarose, cellulose, and derivatives (neutral polymers); and pectin

2. Synthetic polymers, for example, PEG and PLA (polyesters), acrylates, and PVA

3. Combinations of natural and synthetic polymers, for example, alginate and Pluronics^© (block copolymers PPO-PEO-PPO)

Hydrogels can be classified according to several criteria. The most widespread and used one is based on the nature (physical or chemical) of the gel [1, 7, 9]:

1. Physical or “reversible” hydrogels:

   1. Simple entanglement systems, in which the network is held together by molecular entanglements or crystallites

   2. Ion-mediated or “ionotropic” networks, in which the network is stabilized by interaction between polyelectrolyte and multivalent ions of opposite charges (e.g., alginate)

   3. Thermally induced networks, in which the heating (or the cooling) induces the structure formation (e.g., Pluronics^© micelles which require heat to be formed)

1. Chemical or “permanent” hydrogels, mainly covalently bounded

   1. Cross-linked polymers (by radiation, by chemical cross-linkers, and by multifunctional reactive compounds)

   2. Copolymers obtained from monomer + cross-linkers or monomer + multifunctional macromer

   3. Conversion of hydrophobic polymer in a hydrophilic polymer

   4. Polymerization of a monomer in the presence of a solid polymer to obtain an interpenetrating network (IPN)

Other possible classifications are based on the following:

Macromolecular structure [1]:

1. Cross-linked or entangled networks of linear homopolymers, linear copolymers, and block or graft copolymers

2. Polyion–multivalent ion and polyion–polyion or H-bonded complexes

3. Hydrophilic networks stabilized by hydrophobic domains

4. Interpenetrating networks (IPNs) or physical blends

Final form of the system [1]:

1. Solid molded forms (e.g., soft contact lenses)

2. Pressed powder matrices (e.g., tablets for oral delivery)

3. Microparticles (e.g., bioadhesive carriers)

4. Coatings (e.g., on pills and capsules)

5. Membranes (e.g., in transdermal drug delivery patches)

6. Encapsulated solids (e.g., in osmotic pumps)

7. Liquids (e.g., solutions that form gels on heating or on cooling)

1.2. Swelling hydrogels: Nanostructure characterization

The three most important parameters, useful in order to characterize the nanostructure of the hydrogel network, are [2, 3]:

1. The polymer volume fraction in swollen state, ${\varphi }_{2,\text{s}}$

2. The number-average molecular weight of the polymer chain between two neighboring cross-linking points, ${\overline{M}}_C$

3. The corresponding mesh size, $\xi$

In this paper, the notation convention will be 1, solvent (water); 2, polymer; 3, solute (drug); and higher numbers for other components (e.g., ions).

1.2.1. Polymer volume fraction in swollen state

The polymer volume fraction in swollen state, ${\varphi }_{2,s}$, is a measure of the amount of water which can be retained in hydrogel, and it is given by the ratio between the dry polymer volume, $V_P$, and the swollen gel volume, $V_G$. It can be seen also as the reciprocal of the volumetric swelling ratio [13], $Q$ :

${\varphi }_{2,s}=\frac{V_P}{V_G}=\frac{1}{Q}$E1

1.2.2. Molecular weight between two consecutive cross-links

The molecular weight between two consecutive cross-links (of chemical or physical nature), ${\overline{M}}_C$, is an average measure of the degree of polymer cross-linking. In order to estimate its value, the Gibbs free energy of the system has to be evaluated in terms of the elastic and the polymer/water mixing contributions ($\text{Δ}{G}_{total}=\text{Δ}{G}_{elastic}+\text{Δ}{G}_{mixing}$). The equilibrium condition is achieved once the change of chemical potential due to the polymer/water mixing, $\text{Δ}{\mu }_{mixing}$ (expressed on the basis of entropy and heat of mixing), and the change of chemical potential due to the elastic forces, $\text{Δ}{\mu }_{elastic}$, are equal. The change of chemical potential due to the elastic forces can be estimated on the basis of the theory of rubber elasticity [14, 15]. Equating these two terms, an expression for determining the molecular weight between two consecutive cross-links has been derived by Flory and Rehner [16] for a polymer network swollen by absorption of solvent:

$\frac{1}{{\overline{M}}_C}=\frac{2}{{\overline{M}}_n}-\frac{\overline{\upsilon }}{V_1}\frac{\left[ln\left(1-{\varphi }_{2,s}\right)+{\varphi }_{2,s}+{\chi }_{12}{\varphi }_{2,s}^{2 }\right]}{\left[{\varphi }_{2,s}^{1/3 }-{\varphi }_{2,s}/2\right]}$E2

where ${\overline{M}}_n$ is the number-average molecular weight of the polymer, $\overline{\upsilon }$ is the polymer-specific volume, $V_1$ is the molar volume of the solvent, and ${\chi }_{12}$ is the Flory polymer–solvent interaction parameter. For a network obtained starting from a polymer in solution, in a “relaxed” state, with a volume fraction ${\varphi }_{2,r}$, by swelling with further solvent molecules, Flory [17] and Bray and Merrill [18] derived another equation able to estimate the molecular weight between two consecutive cross-links:

$\frac{1}{{\overline{M}}_C}=\frac{2}{{\overline{M}}_n}-\frac{\overline{\upsilon }}{V_1}\frac{\left[ln\left(1-{\varphi }_{2,s}\right)+{\varphi }_{2,s}+{\chi }_{12}{\varphi }_{2,s}^{2 }\right]}{{\varphi }_{2,r}\left[{\left(\frac{{\varphi }_{2,s}^{}}{{\varphi }_{2,r}}\right)}^{1/3}-\left(\frac{{\varphi }_{2,s}}{2{\varphi }_{2,r}}\right)\right]}$E3

The relaxed state is the state of the polymer immediately after the cross-linking but before the swelling. The use of equation 3 has been assessed by Merrill and Peppas [19, 20] working with PVA hydrogels, and it has been used to estimate the interaction parameter, ${\chi }_{12}$, and its dependence upon polymer concentration and temperature [21]. For ionic hydrogels, i.e., hydrogels in which ionic moieties are present, the theoretical treatment has to take into account the Gibbs free energy change due to the ionic nature of the network. Therefore, the change of chemical potential due to the ionic character of the hydrogel has to be taken into account. The problem has been faced out and solved by Brannon-Peppas and Peppas [22], who had derived the following two equations, to describe the swelling of anionic and cationic hydrogels, respectively, prepared in the presence of the solvent:

$\begin{array}{l}\frac{V_1}{4I{M}_r}{\left(\frac{{\varphi }_{2,s}^{}}{\overline{\upsilon }}\right)}^2{\left(\frac{K_a}{{10}^{-pH}+K_a}\right)}^2=={\left[ln\left(1-{\varphi }_{2,s}\right)+{\varphi }_{2,s}+{\chi }_{12}{\varphi }_{2,s}^{2 }\right]}_{1/3}\\ +\left(\frac{V_1}{\overline{\upsilon }{\overline{M}}_C}\right)\left(1-\frac{2{\overline{M}}_C}{{\overline{M}}_n}\right){\varphi }_{2,r}\left[{\left(\frac{{\varphi }_{2,s}^{}}{{\varphi }_{2,r}}\right)}^{1/3}-\left(\frac{{\varphi }_{2,s}}{2{\varphi }_{2,r}}\right)\right]\end{array}$E4

$\begin{array}{l}\frac{V_1}{4I{M}_r}{\left(\frac{{\varphi }_{2,s}^{}}{\overline{\upsilon }}\right)}^2{\left(\frac{K_b}{{10}^{pH-14}+K_b}\right)}^2=={\left[ln\left(1-{\varphi }_{2,s}\right)+{\varphi }_{2,s}+{\chi }_{12}{\varphi }_{2,s}^{2 }\right]}_{1/3}\\ +\left(\frac{V_1}{\overline{\upsilon }{\overline{M}}_C}\right)\left(1-\frac{2{\overline{M}}_C}{{\overline{M}}_n}\right){\varphi }_{2,r}\left[{\left(\frac{{\varphi }_{2,s}^{}}{{\varphi }_{2,r}}\right)}^{1/3}-\left(\frac{{\varphi }_{2,s}}{2{\varphi }_{2,r}}\right)\right]\end{array}$E5

where $I$ is the ionic strength; $K_a$ and $K_b$ are the dissociation constants for the acid and the base, respectively; and $M_r$ is the molecular weight of the repeating unit. The molecular weight between two consecutive cross-links ${\overline{M}}_C$ can be obtained for the two cases solving the suitable equation for the given value of pH.

1.2.3. Mesh size

The mesh size, $\xi$, also known as the correlation length between two cross-links, is a measure of the space available between the macromolecular chains. Therefore, it is related to the ability of a drug molecule to diffuse through the network to be released. It can be calculated by [23]:

$\xi =\alpha {\left({\overline{r}}_0^2\right)}^{1/2}$E6

where $\alpha$ is the elongation ratio (which for isotropically swollen hydrogels is given by $\alpha ={\varphi }_{2,s}^{-1/3}$) and ${\left({\overline{r}}_0^2\right)}^{1/2}$ is the unperturbed root-mean-square of the end-to-end distance for polymer chains between two neighboring cross-links. The last one can be evaluated as follows:

${\left({\overline{r}}_0^2\right)}^{1/2}=l{\left(C_{n}N\right)}^{1/2}=l{\left(C_n\frac{2{\overline{M}}_C}{M_r}\right)}^{1/2}$E7

where $l$ is the length of the bond along the polymer backbone (for vinyl polymers is 0.154 nm), $C_n$ is the Flory characteristic ratio (actually, eq. 7 is the definition of Flory’s characteristic ratio: $C_n={\overline{r}}_0^2/N{l}^2$), and $N=\left(2{\overline{M}}_C/M_r\right)$ is the number of links for the chain. Therefore, the mesh size for isotropically swollen hydrogels can be calculated as follows:

$\xi =l{\varphi }_{2,s}^{-1/3}{\left(C_n\frac{2{\overline{M}}_C}{M_r}\right)}^{1/2}$E8

All these three parameters $\left\{{\varphi }_{2,s},{\overline{M}}_C,\xi \right\}$ are important for the drug delivery processes, since they determine the kinetics of drug diffusion and release. They, which are related to one another, can be determined theoretically or by means of experimental techniques. The ability to tailor the molecular structure of hydrogels (i.e., to design hydrogels with desired values of the three parameters) allows to tailor the hydrogels’ mechanical and diffusive properties.

1.2.4. Structure-properties relationships: Mechanical behavior

From the mechanical point of view, the rubber elasticity theory [24], applicable also to hydrogels for relatively small deformation (less than 20 %), allows to evaluate the tensile stress $\tau$, expressed as the force per unit area of an unstretched sample, swollen by absorption of solvent (analogous to eq. 2):

$\tau =\frac{RT}{\overline{\upsilon }{\overline{M}}_C}\left(\alpha -\frac{1}{{\alpha }^2}\right)\left(1-\frac{2{\overline{M}}_C}{{\overline{M}}_n}\right){\varphi }_{2,s}^{1/3}$E9

For the case in which the polymer is initially in solution at mass concentration ${\rho }_{2,r}$, the tensile stress can be calculated as follows [19]:

$\tau ={\rho }_{2,r}\frac{RT}{{\overline{M}}_C}\left(\alpha -\frac{1}{{\alpha }^2}\right)\left(1-\frac{2{\overline{M}}_C}{{\overline{M}}_n}\right){\left(\frac{{\varphi }_{2,s}}{{\varphi }_{2,r}}\right)}^{1/3}$E10

To characterize the hydrogel, the cross-linking concentration, $c_x$ (sometimes defined as cross-linking density), can be used. The cross-linking concentration is defined as the mole of cross-links for unit volume, and thus, it can be calculated as follows [19]:

$c_x=\frac{1}{\overline{v}{\overline{M}}_C}=\frac{\rho }{{\overline{M}}_C}$E11

where $\rho$ is the polymer density (the inverse of polymer-specific volume). Using the cross-linking concentration and the rubber elasticity theory, the elastic shear modulus at rest, $G_0$, can be evaluated (e.g., equation 4.9c in [15]), in which ${\overline{M}}_C$ can be calculated by eq. 2, 3, 4, or 5:

$G_0=c_{x}RT=\frac{\rho }{{\overline{M}}_C}RT$E12

It is worth noticing that, according to the “equivalent network theory,” in which the network is represented as a collection of spherical “blobs,” whose diameters represent the network mesh size $\xi$ of the entangled structure, it is possible to relate the cross-linking concentration, $c_x$, with the mesh size:

$\xi ={\left(\frac{1}{N_A}\frac{6}{\pi c_x}\right)}^{1/3}={\left(\frac{1}{N_A}\frac{6}{\pi \rho }{\overline{M}}_C\right)}^{1/3}$E13

A comparison between eq. 8 and eq. 13 is not straightforward, since eq. 13 predicts a one-third power dependence of $\xi$ from ${\overline{M}}_C$, whereas eq. 8 requires the knowledge of ${\varphi }_{2,s}$ (plus, if necessary, of ${\varphi }_{2,r}$ and of pH) dependence upon ${\overline{M}}_C$ (given by eq. 2, 3, 4, or 5), and thus, it is not possible to clearly identify the dependence of $\xi$ from ${\overline{M}}_C$. Generally, the equivalent network theory (eq. 13) predicts mesh size smaller than the swelling theory (eq. 8).

1.2.5. Structure-properties relationships: Mass-transport behavior

The release of an active ingredient (AI, a drug molecule) from a hydrogel-based system takes place mainly by diffusion. These systems are prepared by dissolution or dispersion of the drug within the polymer matrix, which is initially dry. Under this condition, the mass transport of the drug through the matrix is negligible; therefore, no release kinetics is observable. When the matrix comes into contact with water (or biological fluids mainly composed of water), the polymer gives origin to a swollen hydrogel network, within which the drug can diffuse and then it can be released (see also section 1.3). The diffusion in hydrogel networks can be concentration gradient-driven (Fickian diffusion) or polymer relaxation-driven (non-Fickian, or viscoelastic, diffusion). The non-Fickian diffusion is known as anomalous transport, and its limiting case has been defined by Alfrey et al. [25] as Case-II transport. The Case-II transport is characterized by a sharp front between the rubbery and the glassy region of the network (see also section 1.3), advancing at a constant velocity, and the rubbery region is a swollen network at equilibrium with the solvent. Because of the constant-velocity advancement of the front, in principle the Case-II diffusion can give rise to a constant mass-flow-rate release of the drug, i.e., a zero-order release kinetics [26].

Vrentas et al. [27, 28] have suggested a simple way to establish the regions in which Fickian and non-Fickian transports take place, on the basis of the diffusional Deborah number, $N_{De.D}$, defined according to eq. 14, in which ${\lambda }_m$ is the characteristic stress-relaxation time of the polymer–solvent system and ${\theta }_D$ is the characteristic time for the diffusion of the solvent in the polymer:

$N_{De.D}=\frac{{\lambda }_m}{{\theta }_D}=\frac{{{\displaystyle \int }}_0^{\infty }sG\left(s\right)ds/{{\displaystyle \int }}_0^{\infty }G\left(s\right)ds}{L_{Ch}^2/D_{1.s}}$E14

The stress-relaxation time, ${\lambda }_m$, can be evaluated by integrals of the shear relaxation modulus, $G\left(t\right)$, over the entire relaxation time spectrum; the diffusion time, ${\theta }_D$, is given by the ratio between the second power of a characteristic diffusional path length (for the solvent), $L_{Ch}$, and the diffusion coefficient of the solvent in the swollen network, $D_{1.s}$ [29]. If the change in solvent concentration during the swelling process is limited, average values for each characteristic time have to be used, and then the full process can be characterized by a single value of the Deborah number. If the change in solvent concentration is large, the Deborah number have to be calculated for both the initial and final stages, and their order of magnitude will be used to characterize the behavior of the system.

Figure 1.

Schematic illustrating the different zones of diffusion, function of the temperature, and penetrant concentration. The solid lines represent lines at constant diffusional Deborah number, $N_{\mathrm{De}.D}$ ; the dashed gray line represents the effective glass transition temperature, $T_g$. The temperatures $T_E$ and $T_V$ represent the temperatures at which the polymer behaves like an elastic solid and viscous fluid, respectively. Figure redrawn on the basis of the suggestions in [27, 29].

The value of the diffusional Deborah number, $N_{De.D}$, discriminates the nature of the diffusive phenomena, as schematized in Figure 1:

1. Large values of the Deborah number ($N_{De.D}\gg 1$) identify Zone I, where the characteristic relaxation time, ${\lambda }_m$, is long with respect to characteristic diffusion time, ${\theta }_D$ : the polymer structure does not change during the water diffusion process, i.e., the polymer remains in its glassy state. The diffusion phenomenon is usually described by the conventional Fick’s law, using coefficients of diffusion constant and independent from water/polymer concentrations.

2. Small values of the Deborah number ($N_{De.D}\ll 1$) identify Zone III, where the relaxation phenomenon is much faster than the diffusion phenomenon. Practically, it is a diffusion through a viscous mixture (the swollen, rubbery hydrogel), a process which can be described by again the conventional Fick’s law, using coefficients of diffusion which are a strong function of water/polymer concentrations.

3. Intermediate values of the Deborah number ($N_{De.D}\approx 1$) identify Zone II, when the two characteristic times are of the same order of magnitude, i.e., the relaxation and the diffusion phenomena take place on the same time scale. This is the transition zone in which the polymer experiences its glass–rubber phase change, the mixture has a viscoelastic nature, and the diffusion is an anomalous transport (and its limit, the Case-II transport), which is non-Fickian.

Several models have been proposed to describe the deviation from the Fickian behavior, among which is the Camera-Roda and Sarti equation [30] that was proven [31] to be able to catch the Case-II and the anomalous transport behaviors with a minimum number of additional parameters. However, despite some experimental evidences of Case II/anomalous transport in polymers used for swellable hydrogel-based delivery systems [32-34], none of the mechanistic models developed to describe drug release have included this feature, considering an instantaneous rearrangement of the polymer chains and therefore limiting their attention to the Zone III (Figure 1) behavior.

Another way to characterize the behavior of the hydrogel network, with particular reference to drug diffusion, requires to consider the so-called swelling interface number, $N_{Sw.I}$, defined [35] by eq. 15, as the ratio between the velocity of the solvent penetration, $v$, and the rate of diffusion, i.e., the ratio between the diffusivity of the solute (the drug) in the swollen network, $D_{3.s}$, and the thickness of the swollen region through which the solute diffusion occurs, $\delta$ :

$N_{Sw.I}=\frac{v}{D_{3.s}/\delta }$E15

The analysis of diffusional Deborah number, $N_{De.D}$, and of swelling interface number, $N_{Sw.I}$, allows to identify the nature of the mass transport. In order to get a zero-order release kinetics, two conditions are expected to be satisfied [29]:

1. The movement of the solvent has to be controlled by polymer relaxation, i.e., $N_{De.D}$ has to be close to unity.

2. The solute diffusion has to be faster than the front movement (solvent movement), i.e., $N_{Sw.I}$ has to be lower than unity.

1.2.6. Modeling of the diffusion behavior

In order to model the diffusive phenomena which take place in hydrogel networks, several models have been proposed and they have been thoroughly reviewed [36-39], and the reader should refer to these reviews in order to have a better view of the problem. Generally, the models useful to predict the diffusivity of a solute “3” in a swollen network “s,” $D_{3.s}$, with respect to the diffusivity of the same solute in the solvent “1,” have the following general form [13] (the ratio $D_{3.s}/D_{3.1}$ is sometimes called “the retardation effect”):

$\frac{D_{3.s}}{D_{3.1}}=f\left(\xi ,{\varphi }_2,r_s\right)$E16

where $\xi$ and ${\varphi }_2$ have been already defined and they are, respectively, the network mesh size and the polymer volume fraction, while the parameter $r_s$ is the size of the diffusing solute. The mechanistic theories which are used in order to build the left-hand side of eq. 16 are known as hydrodynamic theories, obstruction theories, and theories based on free volume. In the following, the basics of each approach are reported along with the most common models.

Hydrodynamic theories

The hydrodynamic theory assumes that the solute molecules, depicted as hard spheres, move through the liquid phase of the network, the diffusion coefficient being dependent upon the drag force exerted by the liquid molecules on the spheres. The binary diffusive coefficient is given by the Stokes–Einstein equation, i.e., $D_{3.1}=k_{B}T/f$, where $f$ is the frictional draw coefficient (for hard spheres of radius $R$ in a liquid of viscosity $\eta$, the frictional draw coefficient is $f=6\pi \eta R$), and the focus of the hydrodynamic theories is on the estimation of the frictional draw coefficient. Cukier [40] proposed eq. 17 for strongly cross-linked networks (rigid polymeric chains, chemical gels) and eq. 18 for weakly cross-linked networks (flexible polymeric chains, physical gels):

$\frac{D_{3.s}}{D_{3.1}}=exp\left[-\left(\frac{3\pi L_c{N}_A}{M_{f}ln\left(L_c/2r_f\right)}\right)r_s{\varphi }_2^{0.5}\right]$E17

$\frac{D_{3.s}}{D_{3.1}}=exp\left[-k_c{r}_s{\varphi }_2^{0.75}\right]$E18

In eq. 17, $L_c$ is the length of the polymeric chain, $N_A$ is the Avogadro number, $M_f$ is the molecular weight of the polymeric chain, and $r_f$ is the polymer fiber radius. In eq. 18, $k_c$ is a parameter depending on the polymer–solvent system. On the basis of the flow through a porous network, Phillips et al. [41] have obtained eq. 19:

$\begin{array}{c}\frac{D_{3.s}}{D_{3.1}}={\left[1+{\left(\frac{r_s^2}{k}\right)}^{0.5}+\frac{1}{3}\frac{r_s^2}{k}\right]}^{-1}\\ k=\left(0.31\right)r_s^2{\varphi }_2^{-1.17}\end{array}$E19

where $k$ is the hydraulic permeability of the medium, considered to be made of straight and rigid fibers oriented in a random three-dimensional way.

Obstruction theories

The obstruction theories are based on the sieve effect due to the presence of an impenetrable polymer network. The diffusion takes place through the holes in the network; therefore, the path length is increased with respect to the diffusion in pure solvent. The retardation effect is thus calculated on the basis of the sieve effect. Different equations are obtained on the basis of the different model for the polymer network. A useful parameter is given by eq. 20, in which $r_f$ is the polymer fiber radius:

$\alpha ={\left(\frac{r_s+r_f}{r_f}\right)}^2{\varphi }_2$E20

Most of the models based on obstruction theories have the same mathematical structure, given by eq. 21. The two parameters $\left\{a,b\right\}$ are listed in Table 1 for different models (for Amsden’s model, $k_1$ is a further parameter depending from the polymer–solvent system).

$\frac{D_{3.s}}{D_{3.1}}=exp\left[-a\cdot {\alpha }^b\right]$E21

Model		$a$	$b$
Ogston et al.	[42]	1.00	0.50
Johansson et al.	[43]	0.84	1.09
Amsden	[38]	$\frac{\pi }{{\left(k_1+2{\varphi }_2^{0.5}\right)}^2}$	1.00

Table 1.

Obstruction theory-based models

By analogy with electrical conduction, Tsai and Strieder [44] proposed a different model, given by eq. 22:

$\frac{D_{3.s}}{D_{3.1}}={\left(1+\frac{2}{3}\alpha \right)}^{-1}$E22

Obstruction theories are better applicable to heterogeneous networks made of rigid polymeric chains.

Combined (hydrodynamic and obstruction) theories

To improve the predictive capability of hydrodynamic and obstruction models, their retardation effects can be multiplied with each other. For example, Johansson et al. [45] combined their obstruction model [43] with the hydrodynamic model by Phillips et al. [41], obtaining eq. 23 which is a better predictor than the two parent models:

$\frac{D_{3.s}}{D_{3.1}}=\frac{exp\left(-0.84{\alpha }^{1.09}\right)}{\left[1+{\left(\frac{r_s^2}{k}\right)}^{0.5}+\frac{1}{3}\frac{r_s^2}{k}\right]}$E23

Another example is given by Clague and Phillips [46], who proposed a structure for the hydrodynamic term, and then they coupled this last with the obstruction model proposed by Tsai and Streider [44]. The resulting model is given here as eq. 24:

$\frac{D_{3.s}}{D_{3.1}}={\left(1+\frac{2}{3}\alpha \right)}^{-1}exp\left[-\pi {\varphi }_2^{\left(0.174\right)ln\left(59.6r_f/r_s\right)\left(r_f/r_s\right)}\right]$E24

Free volume theory

According to the free volume theory, a molecule of solute diffuses through the hydrogel by “jumping” into void which is present in the network. The free volume is the volume of the holes, formed on statistical bases due to random thermal motions. The total free volume is due to the water (solvent, index “ $w$ ” or “ $1$ ”) and to the polymer (index “ $p$ ” or “ $2$ ”), neglecting the presence of the drug, i.e., $v_f={\varphi }_1{v}_{f.w}+{\varphi }_2{v}_{f.p}$. Assuming that the free volume in the polymer, $v_{f.p}$, is negligible, $v_f\approx \left(1-{\varphi }_2\right)v_{f.w}$. The retardation effect can be evaluated on the basis of statistical reasoning, accounting for the probabilities of finding holes in the water free volume, $v_{f.w}$, close enough to the molecule which is moving, accounting that such holes are big enough and accounting for the sieving effects by the polymeric chains. Different estimations for these probabilities give different results for the retardation effect; however, the general structure of models for retardation effect, based on free volume theory, is always given by eq. 25. The two parameters $\left\{a,b\right\}$ are listed in Table 2 for the most common models. The parameter “ $a$ ” plays the role of a sieving factor of the polymer network on the solute molecule, and for drug size smaller than the mesh size $\left(r_s\ll \xi \right)$, $a$ is unity.

$\frac{D_{3.s}}{D_{3.1}}=a\cdot exp\left[-b\left(\frac{{\varphi }_2}{1-{\varphi }_2}\right)\right]$E25

Model		$a$	$b$
Yasuda et al.	[47]	$P_0$	$\frac{B{a}^{*}}{v_{f.w}}$
Peppas and Reinhart	[48]	$k_1\left(\frac{{\overline{M}}_c-{\overline{M}}_c^{*}}{{\overline{M}}_n-{\overline{M}}_c^{*}}\right)$	$k_2{r}_s^2$
Lustig and Peppas	[49]	$\left(1-\frac{r_s}{\xi }\right)$	$\frac{\gamma \pi \lambda }{v_{f.w}}{r}_s^2\cong 1$
Hennink et al.	[50]	$\text{Ψ}$	$k_2{r}_s^2$

Table 2.

Free volume theory-based models

In Table 2, the meanings of the parameters are $P_0$ is the probability of finding an opening between the polymer chains (sieve effect); $a^{*}$ is the effective cross-sectional area of the solute molecule; $B$ is a model parameter; $k_1$ and $k_2$ are two structural constants, different for each polymer–solvent system; ${\overline{M}}_c^{*}$ is a critical molecular weight between cross-links to allow solute passage; $\gamma$ is a numerical factor $\left(0.5\le \gamma \le 1.0\right)$ used to correct overlapping of free volume between more than one solute molecule; $\lambda$ is the jump length, which is roughly equivalent to the solute diameter $\left(\lambda \approx 2r_s\right)$ ; and $\text{Ψ}$ is the sieving factor.

All the models listed in equations 17 to 25 are useful to have an idea of the retardation effect, but still there is not a single approach able to describe the behavior of all the hydrogels. Therefore, for modeling purposes, a phenomenological approach, based on simple equations, with a limited number of fitting parameters, is highly desirable. As a result of the free-volume approach, a simple equation useful to predict the diffusive coefficient for a given molecule specie “i” in hydrogel, with respect to the diffusive coefficient of the same molecule in the swollen network (at equilibrium conditions), $D_{i.H}/D_{i.s}$, has been proposed by Fujita [51], used firstly by Korsmeyer et al. [52], to describe the transport of both the solvent $\left(i=1\right)$ and the solute $\left(i=3\right)$ within the polymeric network:

$\frac{D_{i.H}}{D_{i.s}}=exp\left[-{\beta }_i\left(1-\frac{{\varphi }_1}{{\varphi }_{1.s}}\right)\right]$E26

where ${\beta }_i$ is the single parameter of the equation, and the effect of this parameter on the predictions has been extensively tested by Korsmeyer et al. [53].

1.3. Phenomenology of drug release from swelling hydrogels

In controlled drug delivery, the main goal is to design systems able to give tailored release profiles, as function of the time as well as in response to external stimuli (sensitive systems) and/or in a certain environment (targeted delivery systems). One of the main goals of controlled drug delivery is a system able to give a zero-order release kinetics, i.e., a constant-rate release, for which the mass flow rate of drug released is designed to be equal to the rate of consumption of the drug by the body (by metabolism, excretion, and reactions with its targets), causing a constant plasma concentration of the active principle. Graham and McNeill [54] were the first to prove that by combination of the device geometry, nature of hydrogel, and initial concentration profile of the drug, a zero-order release kinetics is possible working with hydrogel-based systems (morphine release from pessaries shaped as hollow cylinders made of PEO, with a U-shaped initial concentration of the drug through the cylinder wall). This was the starting point for several studies, development, and applications of hydrogel-based pharmaceutical systems for drug delivery. Several pharmaceuticals based on this approach have made their way up to the market.

A full understanding of the phenomena involved, and their detailed mathematical description (modeling), is an objective still unfulfilled. The detailed knowledge of the transport phenomena involved is the key prerequisite in developing a reliable mathematical model useful for the prediction of the release kinetics as function of the formulation parameters or of the external conditions [55]; on the other hand, the availability of a reliable mathematical model of the drug release process could allow to theoretically predict the drug release profile for a given, newly designed, pharmaceutical product, reducing the number of necessary experiments and facilitating the development of new pharmaceuticals [11].

The basics of release mechanisms from swellable hydrogel-based pharmaceuticals are summarized in the following and they are schematized in Figure 2.

Figure 2.

Drug release from a matrix (tablet) made of swellable hydrogel. The phenomena (described in the text) of water uptake, swelling, polymer chain relaxation, drug diffusion, polymer chain disentanglement, and polymer erosion are graphically depicted.

When swellable hydrogel-based matrix comes into contact with water, the following are observed:

1. The water diffuses into the matrix.

2. The water acts as a plasticizer, lowering the polymer glass transition temperature, $T_g$, thus causing the glass–rubber transition, the gel formation, and polymer swelling (polymer chain relaxation).

3. The (soluble) drug dissolves and diffuses through the gel layer.

4. The drug can be released in the dissolution medium.

5. Finally, the polymer dissolves (erodes) at the matrix surface. In case of insoluble drugs, erosion is the main mechanism which allows drug release in the external medium.

This sequence of steps, suggested by Siepmann et al. [56], is coherent with a space distribution of the polymer chains as described by Ju et al. [57, 58], who identified several regions in a hydrogel-based matrix exposed to a water-based dissolution medium:

1. An inner dry glassy core, in which the polymer is practically un-hydrated

2. A glassy gel layer, in which the polymeric chains have very strong entanglements

3. A rubbery gel layer, with highly swollen hydrogels and strong entanglement

4. A diffusion layer, with weak entanglements within the few polymeric chains present

5. The external, or bulk, medium, in which the dissolved polymeric chains are found

These layers, which are not rigidly defined and thus not clearly identifiable, are characterized by a polymer concentration decreasing from the inner regions toward the outer regions, a water concentration decreasing from the outer regions toward the inner regions. Furthermore, they allow to identify several “fronts,” i.e., interfaces which roughly locate the place in which most of the phenomena take place [59]. In particular, the interface between the diffusion layer and the bulk medium is known as the erosion front (because it is the surface on which most of the polymer erosion takes place), the interface between the rubbery gel layer and the diffusion layer is known as the diffusion front (because it is the place in which most of the drug diffusion takes place), and the interface between the glassy gel layer and the rubbery gel layer is known as the swelling front (because it is the place in which the polymer swells, experiencing the glass–rubber transition).

1.4. Swellable hydrogels in drug delivery: A path to understanding

Once the general framework of the phenomena which take place during the drug release from a hydrogel-based pharmaceutical system has been clarified, there are still two needs:

1. The single phenomenon (e.g., the water uptake, the drug diffusion, the polymer dissolution, the shape deformation due to the swelling) has to be quantified, i.e., one or more experimental method(s) have to be pointed out, applied, and validated in order to gather quantitative data useful to confirm the hypothesized phenomenon and to know its extent.

2. Each single observed and quantified phenomenon has to be described by a mathematical sub-model, and all of them together should be collected in a comprehensive mathematical model, able to reproduce quantitatively all the observed phenomena: at this level, the model can be seen as a descriptive tool. Once all the sub-models have been properly tuned using experimental data, the model could become predictive, i.e., able to describe the behavior of a novel system before it is tested or to describe the behavior of an existing system working outside the range of operative conditions within which it has been tested.

The road to the goal of a predictive model is pictorially depicted in Figure 3: on the right of the arrow, several experimental approaches published in literature have been listed (grouped by the main technique adopted in each paper); on the left of the arrow, the succession of modeling approaches proposed in literature is listed (grouped by similarity between the approaches). Of course, the lists are not exhaustive, since a huge number of papers on the topic have been published in the last years. In the remainder of the chapter, the experimental approaches (section 2) and the modeling approaches (section 3) will be reviewed in better details.

Figure 3.

The evolution of experimental approaches (on the right of the arrow) and of theoretical approaches (on the left of the arrow) reported in literature with reference to the study of drug release from matrix made of swellable hydrogels.

2. Experimental approaches

The study of hydrogel-based systems started on simplified systems (simple geometries such as cylinders or sheets, in which only 1D – for example, lateral – water uptake is allowed) and has been continued on more complicated systems (i.e., pharmaceutical forms composed of several units, each with a specific behavior). The complexity of the whole process which takes place during the drug release from hydrogel-based matrices has led to the development of several experimental techniques to deeply understand the behavior of these polymers.

2.1. Factors influencing the drug release

The factors influencing the drug release from hydrogel matrices have been studied for a long time to create a pharmaceutical form with the desired release profile both theoretically and experimentally. In fact, it has to be taken into account that, in addition to the phenomena which take place during a hydrogel hydration, the presence of a local drug concentration influences the polymer behavior [60]. As described in section 1.3, during the water uptake, the polymer chains are subjected to a relaxation, and the presence of an additional component (the drug) changes both the swelling osmotic pressure and the time-dependent relaxation process. Of course, not only the drug has an impact on the release rate from a hydrogel based matrix, but the entire formulation and its constituents influence the matrix behavior.

2.1.1. Influence of the formulation composition

Among these variables, the polymer concentration and its viscosity grade are the most used to control and modulate the drug release [61] due to the fact that the drug release profile is affected both kinetically and mechanistically by these two variables. For example, changes in HPMC (hydroxypropyl methylcellulose) percentages in the formulation can result in a wide range of release rate, but usually, faster drug release rate is obtained with lower HPMC content, as well as faster polymer release is obtained in these conditions. In the analysis of the influence of the HPMC viscosity grade on the drug release, it was demonstrated that the release rate decreases with the increase of the viscosity of the used polymer, but there is a “limiting HPMC viscosity grade” and beyond its value, no further decrease in the release rate is detected; usually, this value is 15000 cP [61]. These analyses can be performed only by evaluating the drug release in the dissolution medium; instead, more complex experimental methods are necessary to evaluate the microscopic variation in surface topography during the swelling such as the cryogenic SEM [62].

In the optimization of the matrix formulation, it is important to evaluate not only the composition of the matrix but even the micrometric properties such as the bulk density of the powders used. The true density of the powders can be measured using a pycnometer. The bulk density can be calculated pouring a pre-weighted amount of powder into a graduated cylinder and measuring the occupied volume. The tapped bulk density can be calculated by measuring the powder volume after tapping the cylinder three times from a defined height and with a constant frequency (i.e., from 2.5 cm height every 2 s); the ultimate tapped density can be calculated after continued tapping till no further volume decrease is observed [62]. The compressibility index was calculated using the bulk and ultimate tapped bulk density.

2.1.2. Influence of the polymer properties

In addition to the type of formulation, the drug release profile from a pharmaceutical form is affected by factors closely related to the polymer, such as its particle size distribution [63]. Studying the release rate of propranolol hydrochloride from HPMC K15M tablet, it was found that the release rate decreased as the particle size of polymer was reduced till about 180 µm. Further reduction caused no further decrease in dissolution rate. The different particle size fractions are usually obtained by sieving. This behavior is probably due to the fact that the larger-sized fraction of HPMC absorbs water faster than the smaller one [63]. These water uptake studies were performed by using DSC (differential scanning calorimetry). Of course, the behavior of the matrix strongly depends on the polymer chosen to constitute the pharmaceutical form. For example, Zuleger et al. [64] found that neither viscosity grade nor the particle size or the compaction pressure of methyl hydroxyethyl cellulose (MHEC) has a strong impact on the release rate of the drug embedded in these matrices. Thus, there is not a general rule from the matrix behavior; it depends on the formulation and drug used. Moreover, even the dissolution medium can influence the release kinetics, both in the case in which the polymer is sensitive to the medium pH and to the ionic strength of the medium.

A recent review on the factors affecting the drug release from hydrophilic matrices was published by Maderuelo et al. [65]. The authors analyzed the effect of different variables, which depend on drug (such as its solubility, particle size, and initial dose), on polymer (such as its particle size, viscosity, and the influence of the dissolution medium), or on formulation (such as manufacturing process) on the drug release kinetics.

2.2. Evaluation of the water uptake, matrix erosion, and drug release

The water uptake is a key parameter to determine the swelling phenomenon of a hydrogel and, as a consequence, the drug release from a matrix. For this reason, several scientists have developed different methods to quantify the water diffusion into the polymeric matrix.

The simplest technique used to evaluate the weight gain after a certain dissolution time, which is directly connected to the water uptake, is the gravimetrical method. According to this method, the weighted pharmaceutical form is immersed in the dissolution medium and, at predetermined time intervals, it is withdrawn from the medium, lightly patted, and weighted, and then it is dried until constant weight is reached [66]. The percent weight gain and the mass loss can be calculated according to the following equations:

$\text{Weight gain }\left(\text{\%}\right)=\frac{\text{wet weight-dry weight}}{\text{dry weight}}\cdot 100$E27

$\text{Mass loss} \left(\text{\%}\right)=\frac{\text{initial weight-dry weight}}{\text{initial weight}}\cdot 100$E28

The percent of mass loss is directly connected to the matrix erosion and drug release (if drug is contained in the polymer matrix) due to the fact that, after drying, only the polymer and drug remain in the matrix, meaning that all the weight loss during the dissolution is imputable to the dissolution of drug and polymer in the external medium. Investigation of matrix hydration and erosion is a fundamental approach to understand the mechanisms involved during the drug release and the parameters which influence the matrix behavior [67]. For example, if the mass loss follows an increasing kinetics with the dissolution time, probably the controlling mechanism in the pharmaceuticals is the erosion of the matrix; on the contrary, if the weight gain continuously increases, the controlling mechanism can be attributed to the swelling. On these bases, Mamani et al. [68] evaluated the influence of components ratio on the release kinetics varying the composition and pH of the dissolution medium.

Thus, the information which can be obtained by the estimation of the mass determination via gravimetrical analysis is the concentration profiles of water and polymer. If the matrix analyzed is composed also of a drug, with the aid of an analytical technique to determine the amount of drug released (i.e., the spectrometric method), all the macroscopic profiles can be measured [69]. Indeed, during the dissolution, the drug leaves the pharmaceutical form and dissolves in the dissolution medium. A typical dissolution test is performed using the standardized apparatuses approved by the USP (United States Pharmacopeia), which are numbered from 1 to 7 [70]. Withdrawing samples of the medium at certain time intervals, it is possible to know the amount of drug released and then, by the difference with the initial drug content, which is usually known by the manufacturer, the amount of drug still contained in the tablet. Repeating these tests for several dissolution times, the drug concentration profile can be evaluated. Usually, the percentage of drug release is evaluated according to the following expression:

$\text{Drug release} \left(\text{\%}\right)=\frac{\text{drug released in the dissolution medium}}{\text{initial drug content}}\cdot 100$E29

The typical drug release profile from a matrix composed of a hydrogel (HPMC) and a drug (theophylline) is shown in Figure 4, red circles on the right axis [34]. As could be expected, the release kinetics is higher at the beginning of the dissolution (after 3 h of dissolution, about 30 % of the drug is released) and tends to decrease with time, till all of the drug contained into the matrix reaches the dissolution medium. On the contrary, the mass of the drug still contained in the matrix is shown in Figure 4 as green stars, on the left axis. It can be measured by the difference between the drug released and the initial amount of drug in the tablet or by a more direct measure. According to this method, after the chosen dissolution time, the tablet is removed from the medium, weighted, and dried. Then, the tablet is immersed into a medium till complete dissolution, and another spectrometric analysis can be performed on the medium to evaluate the concentration of the drug. Usually this technique is used to perform a check on the total amount of drug in the matrix, which can be different with respect to the theoretical one, i.e., due to manufacture errors.

Figure 4.

Mass evolution inside a swelling matrix during a dissolution test. Drug and polymer mass evolutions (in green) are reported on the left axis and water mass evolution (in blue) on the right axis. On the red axis, the percentage of drug released is shown.

Using the gravimetric technique to evaluate the total amount of water absorbed by the matrix, the method previously described to evaluate the weight gained by the pharmaceutical form can be used. After a certain dissolution time, the matrix can be removed from the dissolution medium, weighted, dried, and weighted again. By the difference between the wet and dry weights, the amount of water taken during the dissolution (after the predetermined time) can be evaluated. In Figure 4, the mass of water absorbed by the tablet is shown as blue squares, on the right axis. As could be seen from the graph, the amount of water which can be absorbed by a hydrogel is well beyond the initial amount of polymer, reaching ten times its value during the time interval shown. The last component which has to be quantified is the polymer. This can be done by a simple technique which provides that once the wet weight of the matrix after a certain dissolution time is known and once the drug and water amounts are quantified, the polymer mass can be calculated by the difference between these masses, according to the following expression:

$\text{Polymer mass}=\text{wet weight} \left(\text{after dissolution}\right)- \left(\text{water}+\text{drug}\right)\text{ masses inside the matrix}$E30

Another technique which allows to directly evaluate the amount of polymer dissolved in the dissolution medium and then of the polymer into the tablet is the sulfuric acid technique, recently re-proposed by Ghori et al. [71]. This method consists of taking from the dissolution medium samples of 1 mL at predetermined dissolution times and, after filtration, adding 1 mL of a solution of 5 % phenol in 0.1 M hydrochloric acid, followed by 5 mL of concentrated sulfuric acid. After 10 min of shaking, the solution is cooled (the reaction taking place is strongly exothermic), and a spectrometric analysis can be performed, at 490 nm of wavelength, due to the fact that the solution assumes an intense colored aspect. Besides, this method is more complicated than the measurement of the difference between components; it is a direct measurement of the polymer mass, which means a more reliable matrix characterization.

2.3. Mass evolution inside the matrix

The described gravimetric technique has the practical advantage to be a very simple technique to determine the macroscopic behavior of a pharmaceutical matrix (i.e., the total water uptake, the total mass of polymer eroded by the dissolution medium) or films [72]. On the other hand, it is a destructive technique (different samples have to be removed from the dissolution medium and dried, one for each dissolution time) which gives no information about the concentration profiles inside the matrix and about how they vary with time. For this reason, it is clear that the macroscopic analyses described above are not sufficient to characterize the matrix behavior. Over the years, several experimental techniques have been developed to evaluate the mass evolution inside pharmaceutical matrices.

2.3.1. Image analysis technique

To evaluate the polymer concentration across the gel layer during the tablet hydration, a semiquantitative in situ technique based on the image analysis was developed [73]. According to this method, a HPMC matrix tablet was mounted on a weighted pin and inserted in a vessel filled with the dissolution medium. This stirred vessel was placed on top of a light box. Depending on the orientation of the tablet, the swelling along the radial or axial position can be observed and recorded by a black-and-white camera equipped with a macroscopic zoom lens mounted above the tablet. The light source was represented by a light box containing two fluorescent light tubes. The whole apparatus was closed to avoid the entrance of external sources of light. Each digitalized image had a gray level from 0 (black) to 255 (white) representing the relative scattered light intensity. To understand and explain the results obtained, it was necessary to correlate the scattered light intensity with the concentration of HPMC in the gel layer. This has been done in preparing a series of polymer solutions at known concentrations and evaluating their scattered light intensity, and, by the fitting of these data, a calibration curve to interpret the data was obtained. Due to the complexity involved in the simultaneous evaluation of what happen during the polymer matrix swelling, Bettini et al. [59] used the image analysis to evaluate the movement of the swelling front (as defined in section 1.3) of a tablet during the dissolution in a simplified system which allows the water penetration only in the lateral surface of the tablet. This system was obtained confining the tablet between two transparent Plexiglass disks; in this way, only the radial swelling is allowed, making the system simpler to analyze. Recording the swelling of the matrix during dissolution through the transparent disks, the obtained images were analyzed with image analysis software. Once more, this technique gives no information about the mass fraction evolutions inside the matrix.

On these bases, a nondestructive test to determine the water mass fraction along the radial profile of a pure HPMC tablet was developed by Chirico et al. [55]. Being the tablet composed only by the polymer, this work was focused only on the water penetration inside the matrix, and the tablet was confined between two glass slides to limit the water uptake only in the radial surface. The system composed of the tablet, confined between the glass slides was immersed in the dissolution medium, colored by adding trypan blue to enhance the visibility of the water penetration, and recorded with a camera. Both the camera and the dissolution chamber were closed in a dark room to avoid exposure to external light. Light intensity profiles were evaluated starting from the pictures by the image analysis technique. The analyses were performed, considering the picture as a matrix composed of pixels with intensity values ranging from 0 (white) to 255 (black). To relate the water mass fraction to the light intensity, a relation was adopted:

${\left(\frac{I-I_0}{I_{max}-I_0}\right)}^{\gamma }=\frac{{\omega }_1-{\omega }_{10}}{{\omega }_{1eq}-{\omega }_{10}}$E31

where ${\omega }_{10}$ and ${\omega }_{1eq}$ are the initial water mass fraction in the tablet (in case of a dry tablet, it is equal to 0) and the equilibrium water mass fraction in the hydrogel, respectively. Concerning the light intensity, it was varied between $I_0$, which is a plateau value observed at the matrix core, corresponding to the brightest image, and $I_{max}$, which is the maximum value observed at the fully hydrated gel, corresponding to the darkest image. Since both the variables ($I$ and ${\omega }_1$) vary from 0 to 1, the authors chose to apply the simplest relation between them, a linear one, which is realized when parameter $\gamma$ is equal to 1.

A typical water mass fraction profile obtained via image analysis is reported in Figure 5, a) as gray dotted lines. Moreover, an image of the hydrate tablet after 1 day of dissolution is shown as an inset. Comparing the evolution of the water mass fraction calculated via image analysis and the appearance of the tablet, the correspondence is clear. In fact, in the picture, still visible is the glassy core, in which only a small amount of water is penetrated and the water fraction is the smallest; then it starts to increase due to the fact that the image becomes darker. Then, the swelling front (as defined in section 1.3) becomes visible with a decrease of the light intensity, which leads to an apparently inconsistent decrease of the water mass fraction; actually, this decrease is due to the presence of the swelling front. Then, the light intensity increases gradually, due to the fact that the gel is more hydrated as the radius increases, until the erosion radius is reached, which indicates the final dimension of the swollen tablet, corresponding to a light intensity and water fraction equal to 1. The main drawback of the described technique is the fact that only the water uptake inside the matrix can be evaluated due to the fact that the tablet is composed only of the hydrogel. In fact, the presence of the drug could disturb the image analysis, leading to mistakes in the evaluation of water concentration inside the pharmaceutical form. For this reason, a more complete method taking into account also the drug presence has to be developed.

Figure 5.

a) Water mass fraction along the radial direction after 48 h of dissolution of a radial system. The full squares are the experimental data obtained by the gravimetric analysis; the dotted gray curve is the result of the image analysis; the dashed black curve is the model prediction. In the inset, the real photo of the matrix used for the image analysis. b) Mass fraction evolutions along the radial direction for semioverall system after 6 h of dissolution (the values are averaged along the thickness direction). The scatters are the experimental data; the curves are the model predictions. In both the graphs, the vertical bars represent the standard deviation between the experimental runs, and the horizontal bars represent the extension of the section cut by the punch.

2.3.2. Gravimetric technique

In 2009, Barba et al. [69] proposed a method to quantify all the components inside a matrix at a given dissolution time. This approach, based on the gravimetric technique, was applied firstly on a simplified system, which allows the water uptake only from the lateral surface (the so-called radial system), and then to a more complex system in which the water can penetrate along the whole exchange area (the so-called semioverall system) [34]. The experimental procedure used to evaluate the mass evolutions (polymer, HPMC; drug, theophylline; and water) inside a swollen matrix is illustrated in Figure 6. In the case of the radial systems, the first step was to confine the matrix, containing a polymer, a drug, and eventually an initial water content, between two glass disks. Instead, in the case of the semioverall systems, the matrix was glued at the center of one glass slide, which represents the symmetry plane for the experiment; in fact, the tablet weighted the half part with respect to that used in the radial system. Then, the system was immersed in the dissolution medium until a predetermined time; at the end of this time, the swollen matrix was removed from the dissolution medium (step 2 in Figure 6) and the excess water was removed; in the case of the radial system, the superior slab was carefully removed. The swollen matrix was cut by a series of punches characterized by a thin wall centered in the matrix center (step 3 in Figure 6). The punches had a decreasing diameter: chosen the first punch and centered, all the matrix material external at the punch is recovered on a watch glass (step 4 in Figure 6) and weighted (step 5 in Figure 6). This material contained all the polymer, drug, and water in the wet tablet at greater radii than the punch one. Then, a punch with a small diameter was used for a second cut and, once more, the material recovered in another watch glass. These operations were repeated since a central core was reached, recovered, and weighted. All the samples collected were dried in oven at 45 °C (step 6 in Figure 6) until a constant weight was reached, meaning that all the water initially contained in the wet material was removed (step 7 in Figure 6). These last three phases can be combined together using an infrared moisture meter, as suggested by Tahara et al. [74]. Using this device, the wet material was weighted and dried simultaneously at a temperature of 120 °C, until the ratio of decreasing sample weights becomes less than 0.5 % for 5 min. In the last step of the gravimetrical analysis, each glass was immersed into a vessel containing distilled water till complete dissolution with the aim of quantifying, spectrometrically, the amount of drug still in the matrix.

Figure 6.

Experimental procedure used to evaluate the mass evolutions inside a swollen matrix by gravimetric technique.

Using this method, it is possible to evaluate how the amounts of drug, polymer, and water vary not only with the dissolution time but even with the matrix radius. A typical behavior of a matrix composed initially of 50 % wt of polymer and drug after 6 h of dissolution is reported in Figure 5, b). In this graph, the vertical bars represent the standard deviation between the experimental runs, and the horizontal bars represent the extension of the section cut by the punch. In fact, if the external section is considered, it ranges from 0.6 to 0.8 cm, which means that the punch used had a diameter of 0.6 cm and the wet matrix material recovered had a diameter greater than 0.6 cm and lower than 0.8 cm, which is the maximum radius reached by the swollen matrix. The experimental data were represented at the mean radius of the section of interest, as an average value of mass amount in that section. In Figure 5, the matrix was divided into three sections, two annuli and a central core, and the internal drug and polymer profiles are illustrated as full stars and empty triangles, respectively. As could be seen, both the polymer and drug tend to decrease with the radius; this is easily predictable because, due to the water penetration into the matrix, this tends to unfold the polymeric chains and the drug can easily diffuse outward in the dissolution medium. On the other hand, in the central glassy core, the water is not penetrated yet, and the fractions of the drug and polymer are higher. Concerning the polymer behavior, its fraction decreases with the matrix radius; this is due to the swelling phenomenon which takes place from the interface of the dissolution medium (the erosion front, as defined in section 1.3) toward the inner regions. Concerning the water profile in the swollen matrix, it is shown as full squares in Figure 5, b). The water mass fraction increases with the matrix radius, till it reaches a fraction equal to 1 at the interface with the dissolution medium, which is completely composed of water.

With the aid of this technique, the swollen tablet is completely characterized not only from a macroscopic point of view but also microscopically, because even the internal concentration profiles are known. As previously mentioned, the main drawback of this experimental method is that it is destructive for the sample analyzed.

The water profiles inside the matrix obtained via image analysis are compared with the gravimetric ones in Figure 5, a). To make comparable the results obtained by the gravimetric and the image analysis techniques, the γ value in eq. 31 was optimized to minimize the distance between the two sets of data for a system composed only of polymer and water. The found value is 0.425; thus, it is far from the unity [75], and it depends on the experimental setup. The gravimetric data are the full squares, and the image analysis results are reported as gray dotted line. As could be seen, in the most part of the profile, the data are in good agreement, even if they originate from very different techniques.

2.3.3. Magnetic resonance imaging

Nuclear magnetic resonance imaging (NMRI) is a technique used to produce a two-dimensional map of the density of nuclei inside a thin slice of a more complex object. Besides, the NMRI technique has been extensively used in various fields, such as medical diagnostic and polymer science; in recent years, it has been applied as a noninvasive method to study the behavior of hydrophilic matrices [76]. Usually the nucleus of interest is ¹H in water molecules; thus, the images show the spatial variation of the local water concentration inside the sample; for this reason, this technique is of great aid to characterize the water content inside a swollen matrix during the dissolution. Moreover, the NMRI technique provides an excellent method to evaluate the self-diffusion coefficients, a measure of transport due to the Brownian motion of molecules in the absence of a chemical concentration gradient [77]. The image experiments were usually carried out to measure the water diffusion inside a polymeric matrix and to observe the changing diameter of a tablet during the swelling [78].

Nuclear magnetic resonance imaging (NMRI) is a branch of NMR spectroscopy, usually used to study the swelling and diffusion phenomena of pharmaceutical forms under static conditions (without stirring). The extent of the swelling can be identified by this technique, identifying the position of the interface between the dry core of the tablet and the swollen region. On the other hand, the erosion front can be identified at the interface between the dissolution medium and the swollen region. A great improvement in the use of the MRI technique had been the development of a release cell where the images could be recorded directly at different times during the dissolution, happening in a stirred device [79]. The stirring of the medium was realized in a small release cell, which was equipped with a rotating disk. Gluing the tablet under the rotating disk, the relative movement of the tablet and the dissolution medium ensured the stirring. During the test, the MRI probe was placed at 1 mm below the rotating disk to optimize the tablet position. The release cell is connected to a vessel containing the dissolution medium, kept at constant temperature, which is continuously pumped in the cell by a peristaltic pump. The described setup was used to describe the water uptake into polyethylenoxide (PEO) matrices during a dissolution test [79]. The disk rotation was switched off 2 min before the start of imaging sequences, both in axial and radial directions, and samples of the dissolution medium were collected. The study of the water uptake in the tablet was performed using a multi-spin pulse sequence, in which the signal intensity is weighted on the bases of the proton spin–spin and spin–lattice relaxation times. Due to the very short proton spin–spin relaxation time, differences on the image intensity were mainly dependent on the properties of the dissolution medium. Once the experimental setup was built, it was possible to use this device to better understand the swelling and erosion phenomena happening during the dissolution of a complex matrix. For example, extended release tablets of HPMC were studied via this NMR microimaging technique, and the influence of additives solubility on the matrix behavior was analyzed [80]. By the use of this technique, the authors found that the rate of matrix erosion depends not only on the polymer fraction but even on the solubility of the additives used in the formulation and on the HPMC substituent heterogeneity [81].

In Figure 7, a), a typical example of the use of the NMRI technique to evaluate the water content of a matrix during dissolution is shown. The real picture was taken after 2.7 h of dissolution. At each pixel of the image (128x128) is associated a numerical value that, with the proper scaling factor (given by the instrument), can be related to the proton T2 relaxation. The proton T2 relaxation in turn is mainly attributed to water protons; therefore, it can be related to the amount of water through specific calibration curves. These calibration curves were obtained by analyzing water solutions with known concentrations of solutes. The calibration curve equation obtained with an exponential decay function of the second order was

${\text{H}}_2\text{O} {\text{\%}}_{\left(\text{w}/\text{w}\right)}= a_1-a_2\exp \left(-\frac{\text{T}2}{a_3}\right)-a_4\exp \left(-\frac{\text{T}2}{a_5}\right)$E32

where $a_1÷ a_5$ are specific parameters. The range of confidence of the calibration curve is 60–100 % w/w of water content [34]. The experimental results were shown using different colors to indicate different hydration levels. For example, the blue color was used to indicate the layer of the matrix which showed a percentage of water, included between 60 and 65 %. The results showed a water concentration profile along the “ $r$ ” and “ $z$ ” axes with the 97 % w/w of water content on the erosion boundaries and with a decreasing trend going toward the tablet core. It is worth noting that in this system, the most hydrated layer, 90–97 % w/w, is extremely thin, due to the high system erosion.

Figure 7.

a) On the left, NMR image of a swollen matrix composed of HPMC after 2.7 h of dissolution. At the top, experimental water fraction distribution; at the bottom, calculated water fraction distribution. b) On the left top, photo of a semioverall system after 2 h of dissolution (the matrix was cut along its diameter). On the right top, experimental water mass fraction obtained by texture analysis; on the bottom right, model predictions. In both images, the dimensions are in mm; the color scale is referred to the water content fraction (black = lower water content; gray = fully hydrated matrix).

2.3.4. Texture analysis

A recent approach in pharmaceutical research and development consists of the study of the correlation between drug release and polymer hydration via texture analysis. In general, variation in the textural and physic-mechanical properties could be associated with changes in gel layer and glassy core [82]. One of the first approaches in using the texture analysis to correlate the mechanical properties of the polymeric gels and the water content was that of Yang et al. [82]. They proposed the use of a texture analyzer to evaluate the stress–strain behavior of the polymeric matrix. A matrix composed of HPMC and PEO at different compositions was subjected to dissolution tests. After predetermined time intervals, they were removed from the dissolution medium and subjected to a texture analysis. A probe of 2 mm in diameter was forced to penetrate inside the matrix and a transducer measured the resistance force opposed by the material to the probe penetration. The penetration was ended when a trigger force was reached (0.7 g), the value chosen in a manner in which it was possible to distinguish the glassy polymer from the polymer gel. From the slope of the force–displacement diagrams obtained, it was possible to establish the position of the gel layer; in fact, the resistance encountered by the texture probe is dependent on the gel strength; thus, in the inner part of the tablet, where the water was not penetrated and the gel was strong, the slope of the diagram was very high. On the contrary, in the external parts of the matrix, where the water was penetrated, the resistance encountered by the probe is lower, due to the fact that the gel was hydrated and the slope of the diagrams results to be less pronounced. Based on these considerations, the glassy core and rubbery gel interface can be determined, as well as the overall gel strength. To validate the reliability of these evaluations, the texture experimental data were compared with the position of the gel fronts measured by the NMRI technique. The main drawback of this method is the fact that it is not quantitative: in fact, only the position of the gel front can be quantified, not the water content in the matrix.

In addition to the determination of gel layer thickness and the movement of the erosion and swelling fronts, the texture analysis could be used to determine the total work of probe penetration inside the swollen matrix [66]. Using a method similar to that described previously, swollen thickness was determined by measuring the total probe displacement value (given by the instrument), the axial swelling can be determined by the force–displacement diagrams, and, finally, the work of penetration can be determined by the integration of the force–displacement curve (the area under the curve). The total work of penetration indicates the matrix stiffness or rigidity, and its evaluation can be helpful to understand the effect of an additive on a matrix formulation and how it affects its properties. Moreover, it has to be considered that the gel strength of a matrix system has a great relevance in vivo due to exposure to destructive forces in the gastrointestinal tract.

In 2008, a novel nondestructive probe equipped with a texture analyzer was developed [83], which can be used directly in the dissolution medium during the run. With the use of this probe, the swelling, the erosion, and the positions of the rubbery region and of the glassy core can be simultaneously measured. The probe is constituted by a body through which two supports can move. A circular plate with a central hole contains a fixed probe, which can be exposed outside when the plate is pushed upward. The plate is connected to the movable supports in the probe body. During the analysis, the body of the probe and the plate move downward in the dissolving tablet direction. The loading cell is calibrated in a manner that does not trigger the cell with the movement of the plate in the dissolution medium (the resistance is too low). On the contrary, when the plate reaches the surface of the swollen tablet, it is forced to stop its movement or to move upward if the tablet continues its expansion (due to the swelling). This movement forces the supports connected with the plate to slide into the probe body; this action triggers the load cell which starts to record the resistance force. This movement continues till the probe passes through the hole in the plate and reaches the swollen tablet. The probe is stopped when it reaches a maximum force of 150 g and returns in its original position. This procedure is repeated at preset time intervals and allows the evaluation of the relationship between the thickness of the glassy core and the swollen rubbery region with time. Using the texture analysis technique, various mathematical correlations can be obtained, such as the correlation between the drug dissolution and the polymer hydration obtained by Li and Gu [84]. Modified release tablets were subjected to a conventional dissolution test after their snapping into cylindrical caps, which had a diameter equal to the tablet time. The sample prepared allows the water penetration only along one direction (from the surface exposed to the dissolution medium); in this way, the analysis was significantly simplified. Then, these samples were subjected to a conventional dissolution test and then removed from the dissolution medium and subjected to a texture analysis. The relationship between polymer hydration and drug dissolution was obtained, correlating texture analysis parameters (such as the area under the force–displacement graph), drug release rate, and PEO proportion in the matrix by a mathematical fitting based on both linear and multiple regressions.

With the aim to quantify the amount of water penetrating into a hydrogel-based matrix, using a fast and simple technique, such as the texture analysis, a correlation between the percentage of water content and the penetration work resulting from an indentation test was developed [85]. In this case, matrices composed of HPMC and a model drug were confined between two glass slides, to allow the water uptake only by the radial direction and then immersed in the dissolution medium. After predefined times, the system (composed by the swollen matrix and the slabs) was withdrawn from the medium and the superior slab removed. Due to this sample preparation, the water amount in the matrix changes along the radius (at different radii, there is a different hydration level), but axially, the water content is the same. In particular, at the same radius, the water content is the same in all the matrix thickness. A texture analysis, in particular an indentation test, was performed on these swollen samples. A needle was used to penetrate, at a constant rate, into the matrix, and a diagram force–displacement was obtained. This procedure was repeated at several distances from the center of the matrix, to evaluate how the resistance force changes with the radius of the tablet. For each radius value, a force–displacement diagram was obtained, and, by the calculation of the area under the curve, the work of penetration of the needle to penetrate into the tablet was calculated. In parallel, a gravimetric analysis (described in section 2.3.2) was performed on samples hydrated with the same procedure to obtain the water profiles inside the swollen matrix. The gravimetric data obtained were fitted using a Boltzmann-type equation to obtain a continuous trend of the water mass fraction inside the swollen matrix. Then, the results obtained in terms of water content (evaluated by the gravimetric technique) and the work of penetration (evaluated by the texture analysis) were compared to find a relation between them. A simple relation was obtained:

${\omega }_1\left(W_P\right)=\frac{b_1}{1+{\left(b_2{W}_p\right)}^{b_3}}$E33

where $b_1,b_2, b_3$ are equation parameters, ${\omega }_1$ is the water mass fraction, and $W_p$ is the penetration work. As expected, the penetration work decreases, increasing the water content; this is due to the fact that more hydrated is the matrix, more gel is formed, and less force is necessary to penetrate into the sample. The main drawback of this technique is the fact that this relation between the work and the water content is possible and reliable only if the water content along the matrix thickness is constant (the force–displacement curves are linear); in the cases where the matrix swells both in axial and radial directions, this relation is not so immediate.

Thus, this technique was improved to correlate the water content to the slope of the force–displacement diagrams [86]. This technique was calibrated firstly on the matrix which swells only in the radial direction, since the force increases linearly with the penetration of the matrix. The correlation obtained is

${\omega }_1\left(dF/ds\right)=c_{1}Log\left(\frac{dF}{ds}\right)+c_2$E34

where $c_1$ and $c_2$ are equation parameters, ${\omega }_1$ is the water percentage into the matrix, and $dF/ds$ is the slope of the force–penetration curves. In this way, a calibration relation was obtained and, in each radius of the matrix, it is possible to know the water percentage. At this point, the method was extended to the matrix swollen both in the radial and in the axial direction, the so-called “semioverall” system. Concerning the texture analysis of these systems, the shape of the force–penetration curves was very different with respect to the previous one. In fact, being the water content variable even along the thickness (the water penetration happens both in the radial and axial directions), the measured force increases with the penetration of the needle in a nonlinear behavior. In fact, at the beginning of the test, the needle meets the gel layer formed on the top of the swollen tablet; thus, the force necessary to penetrate this layer is very low. Approaching the core of the tablet, where the water is not penetrated enough, the probe needs a higher force to penetrate the tablet, and the values measured by the texture analyzer increase exponentially. This test was repeated for several radii of the swollen matrix, to evaluate how the force necessary to penetrate varies not only with the thickness but even with the radius of the tablet. This behavior is characteristic of all the curves; the force magnitude varies depending on the radius where the penetration is carried on, since the hydration of the matrix occurs also along the radial direction, thus leaving the matrix center the water content increases; that means that the force necessary to penetrate the matrix decreases. Moreover, it has to be considered that also the thickness of the swollen matrix changes at the changing of the radius position. The experimental data are very well fitted by simple exponential curves. It is thus easy to obtain $dF/ds$ which is, in this case, a function of the penetration distance rather than a single value for each test. Once it is evaluated how the $dF/ds$ varies, for each radius, along the matrix thickness, it is possible to calculate the respective water content using eq. 34.

The measurement of water content in the axial direction can be repeated for each penetration point for all the hydration times considered. Thus, a contour plot can be drawn and compared with a real picture of the semioverall swollen tablet; in Figure 7, b) is shown a matrix after 2 h of dissolution. The picture taken is of a swollen tablet cut along a diameter, in such a way that the glassy core and the swollen gel are clearly distinguishable. The contour plot represents the water content in the tablet evaluated by the described technique, and it is shown in Figure 7, in the first quarter of the picture. Areas containing the same water amount are colored with the same filling, and, in general, the water content decreases while approaching the glass slide of the tablet. It can be noted that there are zones where it is not possible to evaluate the water content, due to both a water content higher than the limit imposed by the obtained relationship applicability (relationship between the force/displacement derivative and the water content) and the hydration in the more external zones, so high that the force opposed by the gel to the penetration is negligible and not easily detectable by the instrument. The legend of the color used and of the water percentage at which they correspond is shown in the third quarter of the picture.

2.4. Characterization of the swelling behavior

In the description of the performances of delivery systems for controlled release, the swelling phenomenon has a central role. For this reason, over the years, several scientists tried to quantify and predict the capability of the matrices to swell and thus release the drug eventually contained inside, usually quantifying the so-called swelling ratio, which corresponds to the weight gain defined in section 2.2. One of the first attempts to describe the swelling phenomenon was that followed by Colombo and his coworkers [87]. The goal of their work was to quantify the matrix relaxation measuring the surface exposed during polymer swelling. To evaluate the influence of the swelling phenomenon on drug release, the amount of drug released was plotted versus the area of the system, at the same time, obtaining a linear relationship. This result indicates a direct dependence between the drug release and the tablet surface. Moreover, the authors proposed a new dimensionless number in order to correlate the matrix expansion during dissolution to the increase of the releasing area. This number is called the swelling area number, $N_{Sw.A}$, and it is defined as

$N_{Sw.A}=\frac{1}{D_{3.s}}\cdot \frac{dA}{dt}$E35

where $dA/dt$ is the rate of the area change during the dissolution and $D_{3.s}$ is the drug diffusion coefficient in the swollen phase. Thus, this number represents the ratio between the matrix surface area expansion rate and the drug diffusion, weighting the contribution of the swelling and the diffusion phenomena on the drug release. The matrix swelling was determined by a measure of the releasing area from pictures of the matrix taken during the dissolution. By this technique, the influence of the presence of a hydrogel matrix coating on the swelling kinetics can be evaluated [88].

Later, the behavior of gel layer thickness in swellable matrices loaded with drugs was studied using a colorimetric technique [89]. This method was based on the dissolution of a radial system for predetermined intervals and the use of the image analysis technique to evaluate the position of the swelling, the diffusion, and the erosion fronts. This method is based on the fact that the used drug was light yellow in the dry state and became orange in solution, with a color intensity dependent on its concentration. These measures are able only to evaluate the radial position of the fronts, being a 1D analysis. Moreover, the water concentration at the glassy–rubber interface ($c^{*}$) (thus at the swelling front, as defined in section 1.3) was calculated starting from the glass transition temperature of the polymer ($T_g$) to the experimental temperature ($T_{exp}$), according to the following expression [29]:

$c^{\text{*}}=\frac{T_g-T_{exp}}{\beta /{\alpha }_f}$E36

where $\beta$ is the contribution of water to the expansion of the polymer and ${\alpha }_f$ is the thermal linear expansion coefficient of the polymer. This equation is usually used in terms of volume fraction of the components because, due to the swelling of the polymer, the volume changes in the swollen tablet. Thus, the water volume fraction can be expressed as water volume per volume of polymer. Thus, from the image analysis, the position of the swelling front can be evaluated and by theoretical consideration, the water concentration at the interface calculated.

As seen in the previous paragraph, the NMRI technique is very useful to monitor the swelling and the erosion phenomena. Proton density- and diffusion-weighted images of the hydrating tablets can be acquired at different intervals, and apparent self-diffusion coefficient maps can be generated from the diffusion-weighted images [90]. The image intensities characteristic of the swollen tablets obtained from the NMRI analysis are usually evaluated by software, i.e., ImageJ.

During the swelling of a hydrogel, the mechanical properties of the gel, and then of the release formulation, change continuously. Thus, one of the analyses which can be performed to evaluate how the swelling influences the matrix behavior is the rheological characterization of the polymer. The effective elastic modulus and gel thickness during the swelling can be evaluated by indentation tests [91]. In a standard indentation test, the needle moves till it enters in contact with the sample and reaches a specified load target. The results in terms of load and displacement are used to compare the effective elastic modulus and gel thickness during the swelling. In particular, Hewlett et al. defined the effective modulus from the slope of the stress–displacement curve at a stress of 10 kPa at which corresponds the effective thickness, that is, the displacement at the point where the effective modulus was calculated. The definition of the effective elastic modulus allows to compare the behavior of different polymers. On the other hand, the viscoelastic materials are usually characterized in terms of dynamic mechanical behavior, when they are subjected to an oscillatory test. The authors performed an oscillatory indentation test to a swollen matrix. The indenter was advanced in contact with the sample. Once a specific load was reached, the indenter velocity was oscillated sinusoidally and the response of the material recorded. Higher loads were used at shorter hydration times and lower loads at longer hydration times because, dependent on the amount of water penetrated inside the matrix, the sample changes its characteristic, becoming harder or softer, respectively. Thus, the oscillatory indentation experiments can be used as a small-scale measure of the rheological behavior of polymers to evaluate the viscoelastic nature of the gel layer as it swells.

In general, for a full description of the swelling phenomenon, a single experimental approach cannot be sufficient to consider all the aspects involved. For this reason, a combined approach of four different methods was proposed [92]. In particular, four different methods can be used: determination of the expansion factor, texture analysis, visual swelling observation, and photomicroscopy. The determination of the expansion factor was performed by positioning the sample into a flat-bottomed test tube and adding a small volume of swelling medium. After the medium was penetrated into the sample, the axial expansion was calculated as the ratio between the measured heights of the swollen and dry tablet. Then, the expansion of the tablets and the development of the gel layer thickness were evaluated using a texture analysis. To evaluate visually the swelling phenomena and thus the water uptake, samples were confined between two disks (as in the radial system) and photographed after predetermined dissolution times using a digital camera to obtain the final dimensions of the tablets. Finally, the swelling of polymer slabs was analyzed by photomicrography with a microscope connected to a digital camera. The combination of all these methods allows a deep analysis of the swelling phenomenon.

As described in this section, a lot of experimental techniques which can be used to clearly understand the swelling kinetics of a polymer and the mobility of the water inside the matrix were developed over the years. Each technique described has its advantages and can be a valid aid for a study purpose. A complete review of all the techniques used to monitor the swelling–erosion behavior was proposed by Huanbutta et al. [93], and the pros and cons of each of them were elucidated.

3. Modeling approaches

3.1. Introduction

The knowledge of the mechanisms that affect drug release, in order to obtain a tailored drug release, has been already pointed out to be essential. A considerable aid in the path of understanding can be given by the mathematical modeling of these physical phenomena, which could highlight the main features of drug release from hydrogel-based delivery systems. The modeling of the physical phenomena involved could help in the development and optimization of these systems, sensibly reducing the time and costs required by trial-and-error procedures. The modeling is rather complex because of the presence of several, synergistic and competing, transport phenomena. Several mathematical models have been proposed in literature, from the simple empirical models to the more sophisticated mechanistic models; the most important are reported in Figure 3 and some of them are discussed in the following.

3.2. Empirical models

The first modeling attempt can be traced back to the semiempirical model of Higuchi [94], where the fractional drug release from an ointment (thin film) was related to the square root of the time:

$\frac{M_t}{M_{\infty }}=k\sqrt{t}$E37

$M_t$ is the cumulative amount of drug released at time $t$. $M_{\infty }$ is the mass of drug released at infinite time (equal to the initial drug loading), and $k$ is a constant reflecting the design variables of the system.

In 1985, Peppas and coworkers [95] proposed a generalization of the Higuchi equation, in which the fractional drug release was related to the nth power of the time. The index “n” is an index function of the drug transport regime and of the shape of the delivery system.

$\frac{M_t}{M_{\infty }}=k{t}^n$E38

For example, for a thin film, the exponent “n” in eq. 38 is equal to 0.5 for Fickian diffusive process (Higuchi’s equation), “n” is equal to 1 for the Case-II transport (swelling-controlled drug release), and it takes intermediate values for intermediate behaviors (anomalous transport).

To distinguish between the relative importance of the Fickian release and the swelling-controlled release, Peppas and Sahlin [96] proposed the model

$\frac{M_t}{M_{\infty }}=k_1{t}^m+k_2{t}^{2m}$E39

where $k_1$, $k_2$, and $m$ are constants. The first term on the RHS represents the Fickian diffusional contribution, $F$, whereas the second term the swelling-controlled (case II) contribution, $R$. The relative importance of the two transport mechanisms can be highlighted by the ratio:

$\frac{R}{F}=\frac{k_2{t}^m}{k_1}$E40

Bettini et al. [88] used this equation to investigate the effect of the HPMC molecular weight on the kinetics of drug release. No significant difference was found within the different grades of HPMC and the value of the ratio $R/F$ was smaller than 0.1, suggesting that the drug transport is mainly Fickian diffusion driven.

Due to their simplicity, these equations have been used countless times to analyze the experimental results, often without taking into account that these equations are empirical in nature, and thus, they are useful only under some restrictions: constant diffusivities and negligible swelling (which is not the case of swellable systems [11, 97]). Thus, these equations should be used as simple fitting tools, giving limited information about the release mechanisms, in particular for complex systems such as the hydrogel-based ones.

3.3. The basics of the mechanistic modeling of the hydrogel-based systems

Let us consider a hydrogel-based system loaded with drug, submerged in a dissolution medium. The hydrogel “H” will be formed by the solvent ($i=1$), the polymer matrix ($i=2$), the drug ($i=3$), and the ions ($i=\mathrm{4...}{N}_C$). To describe the movements of these species along with the hydrogel deformation, proper mass balance equations, coupled with the momentum conservation equation, are needed. Two different approaches can be used to describe this system: the multiphasic model and the multicomponent mixture model. With the first, the hydrogel is seen as constituted by different phases, each characterized by their own mass and momentum conservation equations. With the second approach, instead, the system is considered as made of one phase composed of several components. In Table 3 is depicted a general framework useful for the mechanistic modeling of hydrogel-based systems. Most of the literature works dealing with both the drug release from hydrogels and the hydrogel mechanics could be traced back to this scheme, as shown in the following. In the formulation of the multiphasic approach, it has been considered that the drug is already dissolved. The transient and the inertial terms in the momentum balance equations (Table 3 equations (G-I)) have been neglected.

3.3.1. Free energy of a hydrogel network

In Table 3, eq. (K), it is indicated that the stress tensor and the osmotic pressure of the hydrogel are functions of the Helmholtz free energy that in turn is divided in three terms: elastic, mixing, and ionic free energy. The Helmholtz and the Gibbs free energy for these systems are interchangeable since dG=dA+d(pV) and, at constant pressure, the pressure–volume product does not change significantly in swelling [98].

Multiphasic approach	Multicomponent approach
Mass conservation: $\frac{\partial {\rho }_i}{\partial t}=-\nabla \cdot \left({\rho }_i{v}_i\right)i=\mathrm{1...}{N}_C$ (A)
$$\{\begin{array}{c}\frac{\partial {\varphi }_i}{\partial t}=-\nabla \cdot \left({\varphi }_i{v}_i\right)i=1, 2 \\ \frac{\partial {\varphi }_1{c}_i}{\partial t}=-\nabla \cdot \left({\varphi }_1{c}_i{v}_i\right)i=\mathrm{3...}{N}_C\end{array}$$ (B); ${\varphi }_1+{\varphi }_2\cong 1$ (C)	$\frac{\partial {\rho }_i}{\partial t}=-\nabla \cdot \left({\rho }_i{v}_{\mathrm{mix}}\right)-\nabla \cdot j_{i}i=\mathrm{1...}{N}_C$ (D); or; $\rho \frac{\partial {\omega }_i}{\partial t}=-\rho v_{\mathrm{mix}}\cdot \nabla {\omega }_i-\nabla \cdot j_{i}i=\mathrm{1...}{N}_C$ (E); ${\displaystyle \sum }_{i=1}^{N_C}{\omega }_i=1$ (F)
Momentum conservation: $\nabla \cdot \left({\sigma }_H-p\delta \right)=0$ (G)
$\{\begin{array}{c}0=\nabla \cdot \left(-{\varphi }_{1}p\delta +{\mu }_1\dot{\gamma }\right)+{\varphi }_{1 }RT{\displaystyle \sum }_{i=3}^{N_C}\nabla c_i+f_{12}+{\displaystyle \sum }_{i=3}^{N_C}{f}_{1i} i=1\\ 0=\nabla \cdot \left(-{\varphi }_{2}p\delta +{\sigma }_2\right)-{\varphi }_{1 }\left(z_F{c}_{F}F\right) \nabla \psi +f_{21}i=2\\ 0=-{\varphi }_{1}RT \nabla c_i-{\varphi }_1\left(z_i{c}_{i}F\right) \nabla \psi +f_{i1}i=\mathrm{3...}{N}_C\end{array}$ (H)	$\nabla \cdot \left({\sigma }_{\mathrm{mix}}-p\delta \right)=0$ (I)
Distribution of the ionic species: $${\nabla }^2\psi =-\frac{F}{\epsilon {\epsilon }_0}\left({\displaystyle \sum }_{i=4}^{N_C}{z}_i{c}_i-z_F{c}_f\right), \stackrel{electroneutrality}{\overbrace{{\nabla }^2\psi =0}}$$ (J)
Constitutive equations: $\left\{p_{osm},{\sigma }_H\right\}=g\left(A_{el},A_{mix},A_{ion}\right)$ (K); $p=p\text{'}+p^{osm}$ (L)
$f_{\mathrm{ij}}=-f_{\mathrm{ji}}=f_{\mathrm{ij}}\left(\zeta ,v_i,v_j\right)$ (M)	$\begin{array}{l}{j}_i=-{\rho }_i\left(v_i-v_{\mathrm{mix}}\right)\propto -\\ -\rho D_{i}g\left({\mu }_i, x_i, {\omega }_i, M,p_{osm},\psi ,\mathrm{...}\right)\end{array}$; $D_i=g\left({\rho }_{1}or{\rho }_2\right)$ (N)
Unknowns to solve for:
${\varphi }_1,{\varphi }_2,c_3\mathrm{...}{c}_{N_C}, v_1\mathrm{...}{v}_{N_C}$	${\omega }_i, v_{\mathrm{mix}}$

Table 3.

A general framework for the mechanistic modeling of hydrogel-based systems

*Hydrogel ($H$) = Solvent ($i=1$) + Polymer ($i=2$) + Drug ($i=3$) + Ions ($i=\mathrm{4...}{N}_C$)

The hydrogel can be viewed as formed by an elastomeric network with solvent molecules and ions, and the overall free energy of the system will have to consider the presence and the interactions of these species.

Let us consider a dry elastomeric network in which the stress–strain relation is based on the concept of the finite elasticity theory. This theory relates the deformations of the material to the stresses needed to obtain that deformation through the so-called strain energy density function: $W=W\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)$ where ${\lambda }_i$ represents the stretch ratio in each direction (a material for which the stress–strain relationship derives from a strain energy density function is defined hyperelastic). The strain energy density function corresponds to the Helmholtz free energy per unit volume [99], and it can be derived from statistical theories of the elastomeric network. The aim of statistical models is to provide predictive relationships between the molecular structure and topology of the network and its macroscopic behavior.

Several molecular models for rubber networks have been proposed [98, 100], and the two limiting cases are the “affine” and the “phantom” models. Both these simplified network models give a strain energy density function of neo-Hookean type, i.e., in the following is reported the expression derived from the affine network:

$W^{aff}\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)=A_{el}^{aff}\left({\lambda }_1,{\lambda }_2,{\lambda }_3\right)=\frac{1}{2}{k}_{b}T\left[\frac{{\nu }_{el}}{V_0}\left({\lambda }_1^2+{\lambda }_2^2+{\lambda }_3^2-3\right)-2\frac{{\mu }_{el}}{V_0}\ln \left({\lambda }_1{\lambda }_2{\lambda }_3\right)\right]$E41

where $k_b$ is the Boltzmann constant, $T$ is the temperature, $V_0$ is the volume of the dry polymer, ${\nu }_{el}$ is the number of elastic chains, and ${\mu }_{el}$ is the number of junctions in the network. Equation 41 can be reformulated in terms of Finger deformation tensor $B=F{F}^T$, where $F=dx/dX$ is the deformation gradient (in which $x$ and $X$ are the vectors representing a generic point P in the current and in the reference configuration, respectively):

$A_{el}^{aff}\left(B\right)=\frac{1}{2}{k}_{b}T\left[\frac{{\nu }_{el}}{V_0}\left(I_B-3\right)-\frac{{\mu }_{el}}{V_0}\ln \left(II{I}_B\right)\right]$E42

where $I_B=tr\left(B\right)$ and $II{I}_B=\text{det}\left(B\right)$ are the first and the third invariants of $B$.

Since the hydrogel is characterized by the solvent present, the mixing contribution to the free energy has to be considered. As there is no explicit molecular model for the energy of mixing in gel (cross-linked polymer plus solvent), the functional dependence of the free energy of the mixing is generally assumed to be the same as in a polymer solution (linear polymer plus solvent). The classical treatment is based on the Flory–Huggins model, which is built on a lattice model with the assumption of uniform polymer segment concentration throughout the entire system. The free energy of mixing is given by

$A_{mix}=\frac{RT}{V}\left[n_1\ln \left({\varphi }_1\right)+n_2\ln \left({\varphi }_2\right)+n_1{\varphi }_2{\chi }_{12}\right]$E43

where $R$ is the gas constant, $n_1$ and ${\varphi }_1$ are the number of moles and the volume fraction of the solvent, $n_2$ and ${\varphi }_2$ are the number of moles and the volume fraction of the polymer, and ${\chi }_{12}$ is the Flory–Huggins interaction parameter.

When the network chains contain ionic groups (polyelectrolyte gels), there will be additional forces that affect their swelling properties, like the translational entropy of counterions, Coulomb interactions, and ion pair multiplets. It is evident that the mobile and fixed ions present contribute to the system free energy. Several expressions, based on different assumptions, have been used [101-104]. Quite commonly, the free energy due to ionic charges is divided into mobile ion contribution and fixed charge contribution:

$A_{ion}=A_{ion}^{mob}+A_{ion}^{fix}$E44

Following the Frenkel–Flory–Rehner approach, the free energy of a swelling cross-linked polymer is separable and additive, so that

$A=A_{el}+A_{mix}+A_{ion}$E45

This allows to consider singularly the free energy contributions and to derive constitutive equations for polymeric gels. It is quite common to relate the network stress only to the elastic Helmholtz free energy:

${\sigma }_2=2J^{-1}B\frac{\partial A_{el}\left(B\right)}{\partial B}$E46

where $J=\det \left(F\right)$. The mixing and ionic free energy terms instead are often considered as contributing factors to the osmotic pressure ($p^{osm}$), despite some other approaches that are possible [105].

3.3.2. Multicomponent approach: Mass transfer

In 1999, Siepmann et al. [106] developed a mechanistic model based on purely diffusive mass transport equations (Table 3 equations (D) with $\left({\rho }_i{v}_{\mathrm{mix}}\right)\to 0$ and $j_i=-D_i\nabla {\rho }_i$) to calculate the drug release from HPMC-based matrices. The modeled system was made of three species, water “1”, polymer “2”, and drug “3”, and the domain was the 2D axisymmetric representation of a cylindrical tablet (Figure 8). The water and drug transport equations were solved under the assumptions of:

• No volume contraction upon mixing

• Fast drug dissolution compared to drug diffusion

• Perfect sink condition for the drug

• Strong dependence of the diffusivities (of water and drug) on the hydration level

• Affine deformations

Figure 8.

Schematic of the matrix, in light gray, the mathematical domain on which the transport equations were solved.

The mass transport equations (PDEs: partial differential equations) for water and drug in cylindrical coordinates were numerically (finite difference) solved along with the proper initial and boundary conditions:

$$\{\begin{array}{l}\frac{\partial {\rho }_i}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(D_i\frac{\partial {\rho }_i}{\partial r}\right)+\frac{\partial }{\partial z}\left(D_i\frac{\partial {\rho }_i}{\partial z}\right) i=1,3\\ @t=0,\forall r,\forall z,{\rho }_i={\rho }_{i,0}\\ @r=R\left(t\right),\forall t,\forall z{\rho }_i={\rho }_{i,eq}\\ @z=Z\left(t\right),\forall t,\forall r{\rho }_i={\rho }_{i,eq}\\ @r=0,\forall t,\forall z\frac{\partial {\rho }_i}{\partial r}=0\\ @r=0,\forall t,\forall z\frac{\partial {\rho }_i}{\partial r}=0\end{array}$$E47

where ${\rho }_i$ is the mass concentration of the i^th species, ${\rho }_{i,eq}$ is the interface mass concentration of the i^th species, and $D_i$ is the diffusion coefficient of the i^th species, described by a Fujita-type equation [51] (the general form given by eq. 26):

$D_i=D_{i, eq}exp\left[-{\beta }_i\left(1-\frac{{\rho }_1}{{\rho }_{1,eq}}\right)\right]$E48

where $D_{i, eq}/\exp \left({\beta }_i\right)$ are the values of the effective diffusion coefficients in the dry matrix (${\rho }_1$ = 0) and $D_{i, eq}$ are the values of the effective diffusion coefficients in the fully swollen matrix (${\rho }_1={\rho }_{1,eq}$). The polymer erosion/release was described with the following ordinary differential equation (ODE):

$\{\begin{array}{c}\frac{d{m}_2}{dt}=-k_{er}{A}_{er}\left(t\right)\\ @t=0 m_2=m_{2,0}\end{array}$E49

where $m_2$ and $m_{2,0}$ are the polymer mass and the initial polymer mass, $k_{er}$ is an erosion constant, and $A_{er}$ is the erosion surface, that is, the surface exposed to the external medium. The system deformation was obtained, calculating at each time step the amount of water absorbed in the matrix and the released amount of drug and polymer. Considering affine deformations (the shape is conserved), the swelling and shrinking in the axial and the radial directions were obtained. The model parameters, $D_{i, eq}$, ${\beta }_i$, $k_{er}$ were found from the fitting of a wide set of experimental data. The so-tuned model was able to describe the evolution of the global masses of drug, water, and polymer in the system at several dissolution times.

In 2003, Kiil and Dam-Johansen [107] proposed a 1D model able to describe the radial movements of swelling, diffusion, and erosion fronts along with the hydration and drug release, under the (main) assumptions of drug release only in the radial direction, negligible matrix erosion, and constant diffusivities. Even in this case, the system description was based on the mass transport equation (Table 3 equations (D) with $\left({\rho }_i{v}_{\mathrm{mix}}\right)\to 0$ and $j_i=-D_i\nabla {\rho }_i$) with the addition of physical constraints to model the deformation. In particular, the domain was divided radially in three zones individuated by the front radii: zone 1 (solid part) characterized by the swelling radius, zone 2 (gel plus solid drug) between the diffusion and swelling radii, and zone 3 (gel plus dissolved drug) between the erosion and the diffusion radii. The model was made of one ODE, one algebraic equation, and four PDEs to describe, respectively, the movement of the swelling radius, the erosion front, the diffusion front, the dissolved drug concentration, and the water concentration in zones 3 and 2. Most of the model parameters were taken from literature; others, i.e., diffusivities and dissolution kinetic constants, were adjusted by fitting the experimental data. Despite such a complete approach, the model was not able to describe the swelling that occurs when there is no more solid part (zone 1) and its complexity has downhearted from developing a 2D axisymmetric version.

Inspired by Siepmann’s work, in the 2007, Chirico et al. [55] developed a 1D model to describe the hydration of pure HPMC tablets confined between two glass slabs (Figure 9).

Figure 9.

Experimental tablet dissolution scheme (left) and modeling domain (right)

Considering that the main transport phenomenon is the pseudo-diffusion of water from the medium ($r>R\left(t\right)$) into the tablet through the radial direction (Table 3 equations (E) with $\left(v_{\mathrm{mix}}\right)\to 0$), the water mass balance in cylindrical coordinates was solved numerically (with a finite difference method):

$\{\begin{array}{c}\rho \frac{\partial {\omega }_1}{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left(r\rho D_1\frac{\partial {\omega }_1}{\partial r}\right)\\ @t=0, \forall r, {\omega }_1={\omega }_{1,0}\\ @r=0, \forall t, \frac{\partial {\omega }_1}{\partial r}=0 \\ @r=R\left(t\right), \forall t,{\omega }_1={\omega }_{1,eq}\end{array}$E50

where ${\omega }_1$ is the water mass fraction, ${\omega }_{1,eq}$ is the equilibrium water mass fraction in the gel (0.97 [55]), and ${\text{D}}_1$ is the diffusion coefficient of water in HPMC, modeled with a Fujita-type equation (the general form given by eq. 26) to account for the dependency of diffusivity on the water content:

$D_1=D_{1,eq}exp\left[-{\beta }_1\left(1-\frac{{\omega }_1}{{\omega }_{1,eq}}\right)\right]$E51

where $D_{1,eq}$ is the diffusivity in the fully swollen gel ($@{\omega }_1={\omega }_{1,eq}$), and ${\beta }_1$ is a parameter that tunes the degree of dependency of $D_1$ on the water content.

The system swelling was described, considering the variation of the total mass due to the water inlet and the polymer erosion.

$\begin{array}{c}{m}_1={{\displaystyle \int }}^{\text{}}{\rho }_{\left({\omega }_1\right)}{\omega }_{1}dV\\ \{\begin{array}{c}\frac{d{m}_2}{dt}=-k_{er}{A}_{er}\left(t\right)\\ @t=0 m_2=m_{2,0}\end{array}\\ V\left(t\right)=\pi R^{2}H=\frac{m_1+m_2}{⟨{\rho }_{\left({\omega }_1\right)}⟩}\end{array}$E52

where $m_1$ and $m_2$ are the water and polymer masses in the system at a given time, $k_{er}$ and $A_{er}$ are the erosion constant and the erosion surface (surface exposed to the external medium), $V\left(t\right)$ and $H$ are the system volume and the tablet thickness, and $⟨{\rho }_{\left({\omega }_1\right)}⟩$ is the averaged system density. With this approach, after an initial tuning of the model parameters ($D_{ 1,eq}$, ${\beta }_1$, $k_{er}$), the macroscopic system swelling, the amount of water uptake, and the polymer dissolved were correctly described. In their work, the authors compared, satisfactorily, the model results with an extra experimental data: the water mass fraction profile along the radial direction obtained by analyzing the normalized light intensity of swollen tablet pictures; an example of this model application is shown in Figure 5a. The model abilities were later on improved [108] to describe the drug release kinetics from different shaped matrices, bringing back everything to a 1D model valid for slabs, spheres, infinite cylinders, and finite cylinders (where the approach proposed by Coviello et al. [109] was used).

Despite that these transport models were validated against a wide set of experimental data, they were limited by the strong assumption of affine deformation that impedes the use of such models in more geometrically complex systems where the non-affine deformation can sensibly affect the drug release.

In 2011, Lamberti et al. [110] proposed a 2D axisymmetric model, in some ways similar to the Siepmann’s model but with a big step forward: the assumption of affine deformation was no longer needed.

The transport equations for water and drug were solved for the respective mass fraction with a finite element method (FEM), along with the proper initial and boundary conditions (Table 3 equations (E) with $\left(v_{\mathrm{mix}}\right)\to 0$):

$$\{\begin{array}{c}\rho \frac{\partial {\omega }_i}{\partial t}=\nabla \cdot \left(\rho D_i\nabla {\omega }_i\right) i=1,3\\ @t=0,\forall r,\forall z,{\omega }_i={\omega }_{i,0}\\ @r=R\left(t\right),\forall t,\forall z {\omega }_i={\omega }_{i,eq}\\ @z=Z\left(t\right), \forall t,\forall r {\omega }_i={\omega }_{i,eq}\\ @r=0, \forall t,\forall z \frac{\partial {\omega }_i}{\partial r}=0\\ @z=0, \forall t,\forall r \frac{\partial {\omega }_i}{\partial z}=0\end{array}$$E53

where ${\omega }_i$ is the mass fraction of the i^th species, ${\omega }_{i,eq}$ are the interface mass fraction of the i^th species, and $D_i$ is the diffusion coefficient of the i^th species, described by a Fujita-type equation (see eqs. 48 and 51). The constitutive equation for the system density was obtained considering the summability of the specific volume of the single species:

$\frac{1}{\rho }=\frac{{\omega }_1}{{\rho }_{10}}+\frac{(1-{\omega }_1-{\omega }_3)}{{\rho }_{20}}+\frac{{\omega }_3}{{\rho }_{30}}$E54

where ${\rho }_{10},{\rho }_{20},{\rho }_{30}$ are the pure species densities.

The peculiarity of this approach was the decomposition of the inlet water flux to describe the system swelling. In particular, it was considered that part of the total inlet water mass flux was responsible for the tablet swelling, $j_{1,\mathrm{swe}}$, whereas the rest ($j_{1,\mathrm{diff}}$) was responsible for the inner layer hydration. Assuming that

$$\begin{array}{c}{j}_{1,\mathrm{swe}}=k_{swe}{j}_{1,\mathrm{diff}}\\ \text{and}\\ j_1=j_{1,\mathrm{diff}}+j_{1,\mathrm{swe}}=\left(1+k_{swe}\right)\rho D_1\nabla {\omega }_1\end{array}$$E55

where $k_{swe}$ is a new fitting parameter.

Figure 10.

Representation of the water mass balance on a boundary element

The rate of deformation was obtained performing a water mass balance on a boundary element (Figure 10):

$\begin{array}{c}\rho \frac{d}{dt}\left(A\delta {\omega }_{1,eq}\right)=A\cdot j_{1,\mathrm{diff}}-A\cdot j_1\\ v_{\mathrm{swe}}=\frac{d\delta }{dt}=-\frac{j_{1,\mathrm{swe}}}{\rho {\omega }_{1,eq}}=-\frac{k_{swe}{j}_{1,\mathrm{diff}}}{\rho {\omega }_{1,eq}}\end{array}$E56

where $v_{\mathrm{swe}}$ is the swelling velocity of the boundary element. On the other hand, the polymer erosion could lead to a volume decrease; therefore, an opposite constant velocity (the erosion velocity $v_{\mathrm{eros}}=-k_{eros}$) was considered, so that the effective velocity of the boundary element was given by

$v=v_{\mathrm{swe}}+v_{\mathrm{eros}}$E57

The application of this boundary velocity through an ALE (Arbitrary Lagrangian–Eulerian) moving mesh method was proved to describe the matrix non-affine swelling. The model, after an initial tuning of the parameters, was successfully compared to macroscopic results in terms of overall water, drug, and polymer masses, erosion radius, and semi-thickness as well as a picture of swollen tablets.

The approach used in this model overcame the affine deformation at the cost of an additional parameter: the $k_{swe}$ (1° drawback). Moreover, the mass fraction constraint (Table 3 equation (F)) was not implemented along with the PDEs (as well as in Siepmann’s model), possibly leading to unrealistic results in some domain points (2° drawback).

Kaunisto et al. developed a mechanistic model (in 2010 for the dissolution of pure poly(ethylene oxide) tablets [111], extended in 2013 [112] to HPMC matrices loaded with a poorly soluble drug) to analyze the drug release behavior in swellable systems. In this model (considered at constant density), the polymer mass fraction description was coupled with the transport equations for drug and water through the mass fraction constraint (Table 3 equations (E-F) with $\left(v_{\mathrm{mix}}\right)\to 0$), therefore overcoming the Lamberti model 2° drawback. The transport equations were based on a simplified version of the generalized Fick equation [113]. Despite the elegant approach, all the multicomponent interactions, except those with the solvent, were assumed to be zero, and the multicomponent Fick diffusivities were interpreted as “pseudo-binary.” For the water–polymer diffusivity, a “Fujita-type” form was used (eqs. 48 and 51). Even in this case, like in Lamberti’s model, the swelling was described through an ALE moving mesh method, but the swelling velocity was derived from a polymer/solid drug mass balance on the erosion boundaries, with zero additional parameters, therefore overcoming also the 1° Lamberti model drawback.

Recently, in 2015, Caccavo et al. [114, 115] developed a mechanistic model based on the diffusion of water and drug in a concentrate system, accounting for the molar mass gradient, coupled with an ALE method.

Figure 11.

Computational domain for the model of Caccavo et al. [114, 115]

Under the assumption of:

• No volume change upon mixing.

• The drug dissolution within the matrix is fast compared to drug diffusion.

• Perfect sink conditions for drug and constant critical solvent concentration on the erosion front are maintained.

• The water convection contribution to mass transfer is negligible.

• The swelling is due to the water uptake and to the translocation of polymer. The polymer flow at the interface between the tablet and the external medium causes the swelling and the shape change, without polymer release.

• The erosion is due to external fluid dynamics and thus to interactions between the tablet surface and external fluid, which remains constant during the dissolution runs.

The water and drug balances (transport equations) (Table 3 equations (E-F) with $\left(v_{\mathrm{mix}}\right)\to 0$) in cylindrical coordinates (Figure 11) were:

$$\{\begin{array}{l}\rho \frac{\partial {\omega }_i}{\partial t}= \nabla \cdot \left(\rho D_i \nabla {\omega }_i+\rho \frac{{\omega }_i}{M}{D}_i\nabla M\right)\\ {\omega }_2= 1-\left({\omega }_1+{\omega }_3\right)\\ @ t=0 \forall x ϵ\mathrm{\Omega }{\omega }_i={\omega }_{i0}\\ @ xϵ\mathrm{\Gamma }1 \forall t>0 j_i=0\\ @ xϵ\mathrm{\Gamma }2 \forall t>0 j_i=0\\ @ xϵ\mathrm{\Gamma }3 \forall t>0{\omega }_i={\omega }_{i,eq}\\ @ xϵ\mathrm{\Gamma }4 \forall t>0{\omega }_i={\omega }_{i,eq}\end{array}$$E58

where $M$ is the average molar mass. The ideal mixing rule (like in Lamberti’s model) was used for the system density, and the diffusion coefficients were described with Fujita-type equations (eqs. 48 and 51). The domain deformation was described with the ALE method imposing the swelling and the erosion velocities on the erosion boundaries:

$\{\begin{array}{l}\frac{{\partial }^2}{\partial R^2} \left(\frac{\partial r}{\partial t}\right)+\frac{{\partial }^2}{\partial Z^2} \left(\frac{\partial r}{\partial t}\right)=0\\ \frac{{\partial }^2}{\partial R^2}\left(\frac{\partial z}{\partial t}\right)+\frac{{\partial }^2}{\partial Z^2}\left(\frac{\partial z}{\partial t}\right)=0\\ @ t = 0 r=r_0; z=z_0\\ @x ϵ \mathrm{\Gamma }1 \forall t > 0 dr=0\\ @x ϵ \mathrm{\Gamma }2 \forall t > 0 dz=0\\ @ x ϵ \mathrm{\Gamma }3 \forall t >0 \frac{\partial z}{\partial t}=\left(v_{\mathrm{swe}}+v_{\mathrm{er}}\right){\text{}|}_z\\ @x ϵ \mathrm{\Gamma }4 \forall t >0 \frac{\partial r}{\partial t}=(v_{\mathrm{swe}}+v_{\mathrm{er}}{\text{})\text{}|}_r\end{array}$E59

where the uppercase and the lowercase indicate the domain coordinates in the reference and current configuration, respectively. The swelling velocity was derived through a polymer mass balance on the erosion front [114]. The polymer flow toward the erosion front was intended to be used to form new layers of gel, contributing in this way to the swelling phenomena.

$v_{\mathrm{swe}}=-\frac{\left(j_1+j_3\right)}{\rho {\omega }_2}$E60

whereas the erosion velocity was accounted using a constant value $\left|v_{\mathrm{er}}\right|=-k_{er}$.

Thanks to this approach, the evolution of the water, drug, and polymer masses was described, along with the mass fraction profiles inside the swollen tablets (Figure 5b) and the tablet shape evolution (Figure 12).

Figure 12.

Calculated hydration level, in terms of water mass fraction, and swelling at several dissolution times. In the pictures is represented the half portion of the tablet.

Moreover, the model was successfully compared to experimental data of water mass fraction distribution obtained with the texture analysis (Figure 7b) [115]. A modified version of the model [34], with the inclusion of an excipient (lactose), was compared to hydration and swelling data obtained from an MRI technique, proving that in case of Fickian transport, the model is still reliable (Figure 7a) with some discrepancy at higher dissolution times where the fluid dynamic forces, in the specific testing tool, take over and deform the weak gel.

3.3.3. Multicomponent approach: Mass and momentum transfer

Achilleos et al. [101, 116] proposed a modified transport model for a polyelectrolyte gel made of four components (water “1”, polymer “2”, Na⁺ “4”, Cl^- “5”), where the volume transition of the hydrogel was still considered to be driven by the ionic diffusion, but the osmotic pressure is contributed by the mixing, elastic, and electrostatic free energies. Since the system is considered monophasic, with the presence of several species, this kind of formulation can be considered, according to the scheme of Table 3, a multicomponent approach.

Following the mixture theory approach, the mass and momentum balance equations can be written for the mixture, considering a mass average velocity ($v_{\mathrm{mix}}$). The system is considered in isothermal conditions so that there is no need for an energy balance. Along with the momentum and mass balance for the mixture, the species mass balance has to be satisfied. All these equations were nondimensionalized and written in terms of polymer velocity (and not the mass average velocity), since it is the one that enters in the stress constitutive expression:

$$\{\begin{array}{l}{\rho }^{\text{*}}\frac{D^{\left(2\right)}{\omega }_i}{Dt}-{\rho }^{\text{*}}\frac{{\omega }_i}{{\omega }_2}\frac{D^{\left(2\right)}{\omega }_2}{Dt}+\nabla \cdot \left(j_i^{\text{*}}-\frac{{\omega }_i}{{\omega }_2}{j}_2^{*}\right)=0, i\ne 2\\ \nabla \cdot v_2^{\text{*}}+\frac{1}{{\rho }^{\text{*}}}\frac{D^{\left(2\right)}{\rho }^{\text{*}}}{Dt}+\frac{1}{{\omega }_2}\frac{D^{\left(2\right)}{\omega }_2}{Dt}=0\\ \nabla \cdot {\sigma }_2^{\text{*}}-\nabla {p^{\prime }}^{\text{*}}=0\\ {\nabla }^2{\psi }^{\text{*}}\propto -\left({\displaystyle \sum }_{i=4}^{N_C}\frac{z_i{\omega }_i{M}_1}{M_i}-\frac{z_2{\omega }_{2}a{M}_1}{M_{2,m}}\right)=0\end{array}$$E61

where ${\rho }^{*}$, $v_2^{*}$, ${\text{p}}^{\text{'*}}$, and ${\psi }^{*}$ are the dimensionless mixture density, polymer velocity, hydrostatic pressure, and electrical potential, respectively. ${\omega }_i$, $j_i^{*}$, $M_i$, and $z_i$ are the weight fraction, the dimensionless diffusive flux, the molecular weight, and the charge of the i^th species, respectively. The variable $a$ and $M_{2,m}$ represent the polymer network ionization ($0\le a\le 1$) and the molecular weight of a monomer. $D^{\left(2\right)}/Dt$ is the substantial derivative following the polymer velocity.

Equation 61 represents (in order) the species mass balances ($N_C-1$) (Table 3 eq. (E)), the continuity equation (Table 3 sum on $N_C$ of eq. (D)), the momentum balance (Table 3 eq. (I)), and the electroneutrality condition (Table 3 eq. (J)). These equations solve for ${\omega }_i$, ${\text{p}}^{\text{'*}}$, $v_2^{*}$, ${\psi }^{*}$, respectively. It is worth to specify that the mixture stress tensor, with this approach, is equal to the polymer network stress tensor: ${\sigma }_{mix}={\sigma }_2$.

3.3.3.1. Constitutive equations

The diffusive flux of the i^th components is given by (Table 3 eq. (N)):

$j_i^{*}=-D_i^{\text{*}}{\omega }_i\left(\nabla p_i^{osm}{}^{\text{*}}+\frac{M_i}{M_1}\nabla \cdot {\sigma }_2^{\text{*}}+z_i\nabla {\psi }^{\text{*}}\right)$E62

where $p_i^{osm}{}^{*}$ and $D_i^{*}$ are the dimensionless osmotic pressure and diffusion coefficient ($D_1$ was considered function of the water concentration) of the i^th components. The Helmholtz free energy of the mixture, consisting of mixing ($A_{mix}$), elastic ($A_{el}$), and electrostatic ($A_{ion}$) contributions, can be used to obtain the osmotic pressure and the network stress (Table 3 eq. (K)). The mixture and the electrostatic contribution were accounted in the osmotic pressure:

$\begin{array}{l}{p}_i^{osm}{}^{\text{*}}=\left(ln({\omega }_i)+{\chi }_{i,2}{\omega }_2+1\right)-{\omega }_2{\displaystyle \sum }_{\begin{array}{c}\beta =1\\ \beta \ne 2\end{array}}^{N_C}\frac{{\chi }_{\beta ,2}{\omega }_{\beta }{M}_i}{M_{\beta }}-{\displaystyle \sum }_{\begin{array}{c}\beta =1\\ \beta \ne 2\end{array}}^{N_C}\frac{{\omega }_{\beta }{M}_i}{M_{\beta }}\\ +C_{DH}\left(I_s^{1.5}\frac{M_i}{M_1}-3 z_i{I}_s^{0.5}\right)+C_{mf}{a}^2{\omega }_2\left(-\frac{z_i^2{M}_1}{I_s{M}_i}+1\right)\end{array}$E63

where ${\chi }_{i,2}$ is the Flory–Huggins interaction parameter and $C_{DH}$ and $C_{mf}$ are two determined constants. The osmotic pressure is therefore composed of a solvent–polymer mixing contribution, a mobile ions contribution, and fixed charges contribution. The elastic contribution of the Helmholtz free energy was used, according to eq. 46, to derive the network stress tensor:

${\sigma }_2=2 J^{-1}B\frac{\partial A_{el}\left(B\right)}{\partial B}=G_0\frac{{\varphi }_2}{{\varphi }_{20 }}\left(B-\delta \right)=G\left(B-\delta \right)$E64

where $G_0$ is the hydrogel modulus and ${\varphi }_2$ and ${\varphi }_{20 }$ are the volume fraction of the polymer and the volume fraction of polymer at the gelation (stress-free state) so that the shear modulus is $G=G_0{\varphi }_2/{\varphi }_{20}$. The network stress was related to the polymer velocity $v_2$ :

$\frac{D^{\left(2\right)}{\sigma }_2}{Dt}-{\sigma }_2\cdot \nabla v_2-\nabla v_2^T\cdot {\sigma }_2+{\sigma }_2\nabla \cdot v_2=G_0\frac{{\varphi }_2}{{\varphi }_{20 }}\left(\nabla v_2+\nabla v_2^T\right)$E65

These last equations, coupled with the momentum balance, allow to solve for ${\sigma }_2$.

Applying the proper initial and boundary conditions, this model solution (obtained with a FEM method) provided the hydrogel water concentration, deformation, stress, and potential fields as functions of external and gel parameters: solution molarity, degree of polymerization, shear modulus, and problem geometry.

3.3.4. Multicomponent approach: Monophasic theory

On the basis of the poroelastic theory, Suo and his group developed a framework (the monophasic theory) for modeling the mass transport of small molecular species in a three-dimensional polymeric network coupled with large deformations [105, 117]. Once again, the model can be classified, according to the scheme of Table 3, as a multicomponent approach. Indeed, the system is monophasic and the solvent transport is modeled with a diffusive mass transport equation.

The basic idea is that because of its microstructure, a gel is capable of two basic modes of deformation. First, the stretching of polymer chains accompanied by local rearrangement of small molecules allows a gel to change its shape rapidly but maintain a constant volume. The gel behaves just like an incompressible elastomer in this mode, showing a viscoelastic nature. The second mode, which involves long-range transportation of the small molecules, is usually slow and size dependent. The deformation in the second mode allows the gel to swell or shrink in volume, emphasizing the poroelastic nature.

The first mode of deformation was neglected in the earliest publication [117], considering an instantaneous rearrangement of the polymer chain, therefore not considering the viscoelastic nature of the hydrogel. In the wake of this framework, several other works have been published [105, 118-122]. Later on, in seek of completeness, the viscoelastic nature of hydrogel has been then added to the framework [123].

In the following the viscoporoelastic model [123] will be reported, using the “index notation” and the summation convention (with the uppercase and with the lowercase are denoted the variables in the reference and current state, respectively). Let $\mathrm{\Omega }$ be the volume of each solvent molecule in the gel. Taking the dry state of the polymer network to be the reference, it is possible to relate the nominal solvent concentration $C$ to the volumetric expansion $\text{det}F$ as

$\text{det}F:=\frac{V_{tot}}{V_{dry}}=1+\frac{V_{solv}}{V_{dry}}=1+\frac{V_{solv}}{N_{solvmolecules}}\frac{N_{solvmolecules}}{V_{dry}}=1+\mathrm{\Omega }C$E66

As usual, the gel Helmholtz free energy is considered as the sum of the elastic (affine network) contribution and the polymer–solvent mixing contribution. A peculiarity of this approach is the reformulation of the Flory–Huggins $A_{mix}\left(C\right)$ as $A_{mix}\left(F\right)$, using eq. 66, therefore considering the strain energy density function given by the sum of the elastic and mixing free energy.

3.3.4.1. Nonequilibrium thermodynamics

The second law of thermodynamics dictates that the total free energy of a closed system keeps constant or decreases in any physical process. Assuming that the energy is dissipated mainly through two processes, solvent migration $\dot{W_s}$ and viscous deformation $\dot{W_V}$ ($\dot{W_s}$ being dependent on the true solvent flux $j_1$ and $\dot{W_V}$ being dependent on the deviatoric part of the strain rate tensor), it is possible to write

$\frac{\partial \mathrm{\Phi }}{\partial t}-\frac{\partial W}{\partial t}+\underset{V}{{\displaystyle \int }}{\dot{W}}_{s}dV+\underset{V}{{\displaystyle \int }}{\dot{W}}_{V}dV=0$E67

where the first term is the rate of variation of the total free energy ($\mathrm{\Phi }$) of the hydrogel $\partial \mathrm{\Phi }/\partial t=\underset{V}{{\displaystyle \int }}\partial A/\partial t dV$. The second term is the total work ($W$) done to the hydrogel, including the mechanical work by the surface traction and the chemical work via absorption of solvent molecules $\partial W/\partial t=\underset{\partial V}{{\displaystyle \int }}{t}_i\partial x/\partial t dS+\underset{\partial V}{{\displaystyle \int }}\overline{\mu } IdS$, where $t_i$ is the stress acting on the gel surface, $\overline{\mu }$ is the chemical potential of the solvent in the external medium, and $I$ is the molecular flux. The third and fourth terms represent the free energy dissipation due to solvent migration ($\dot{W_s}$) and viscous deformation ($\dot{W_V}$).

Moreover, the system has to respect the solvent molecule conservation $\partial C/\partial t=-\left(\partial J_K\right)/\left(\partial X_K \right)$ that rewritten in terms of $C=C\left(\text{F}\right)$ can be added (using a Lagrange multiplier ($\mu$)) to eq. 67:

$\begin{array}{c}\underset{V}{{\displaystyle \int }}\frac{\partial A}{\partial F_{iK}}\frac{\partial v_i}{\partial X_K}dV-\underset{\partial V}{{\displaystyle \int }}{t}_i\frac{\partial x}{\partial t}dS-\underset{\partial V}{{\displaystyle \int }}\overline{\mu }IdS+\underset{V}{{\displaystyle \int }}{\dot{W}}_{s}dV+\underset{V}{{\displaystyle \int }}{\dot{W}}_{V}dV+\\ +\underset{V}{{\displaystyle \int }}\mu \left[\frac{detF}{\mathrm{\Omega }}{H}_{Ki}\frac{\partial v_i}{\partial X_K}+\frac{\partial J_L}{\partial X_L}\right]dV=0\end{array}$E68

where the tensor $H=F^{-1}$. With the aid of integration by parts and the divergence theorem, it can be shown that to satisfy eq. 68,

$\{\begin{array}{c}\frac{\partial P_{iK}}{\partial X_K}=0\\ P_{iK}{N}_K=t_i\\ \frac{\partial C}{\partial t}+\frac{\partial J_L}{\partial X_L}=\frac{detF}{\mathrm{\Omega }}{H}_{Ki}\frac{\partial v_i}{\partial X_K}+\frac{\partial J_L}{\partial X_L}=0\\ J_K=-\frac{D}{k_{b}T}\frac{detF-1 }{\mathrm{\Omega }detF}{H}_{Ki}{H}_{Li}\frac{\partial \mu }{\partial X_L}\\ \mu =\overline{\mu }\end{array}$E69

where $P_{iK}$ is the nominal stress (first Piola–Kirchhoff stress that relates to the Cauchy stress by $${\sigma }_{\mathrm{mix}}=P{F}^T$$.) defined as

$\begin{array}{l}{P}_{iK}=N{k}_{b}T\left(F_{iK}-H_{iK}\right)+\frac{k_{b}T}{\mathrm{\Omega }}\left[\ln \left(1-\frac{1}{detF}\right)detF+1+\frac{\chi }{detF}\right]H_{Ki}\\ +\eta detF H_{Lj}\left(\frac{1}{2}{H}_{Kj}\frac{\partial v_i}{\partial X_L}+\frac{1}{2}{H}_{Ki}\frac{\partial v_i}{\partial X_L}-\frac{1}{3}{H}_{Ki}\frac{\partial v_i}{\partial X_L}\right)-\mu \frac{detF}{\mathrm{\Omega }}{H}_{Ki}\end{array}$E70

Equation 70 represents the momentum (Table 3 eq. (I) in terms of nominal stress) and the solvent mass conservation equations (Table 3 eq. (D) in the reference configuration) in the hydrogel, derived from nonequilibrium thermodynamics considerations, along with their natural boundary conditions rising from the transformation from the strong to the weak form. The expression of the stress tensor in eq. 70 results in a more complex equation, i.e., with respect to Achilleos’s model, since this accounts for the network elasticity, the polymer–solvent mixing, and the energy dissipations.

3.3.5. Multiphasic models

With the “multiphasic models,” the hydrogel is divided into two (biphasic) or three (triphasic) phases, including a solid polymer network phase, a fluidic solvent phase, and a mobile ion phase for polyelectrolyte hydrogels. In the following will be analyzed the biphasic (polymer plus water) model of Birgersson and Kurnia [124-126]. It has to be said that with the multiphasic approach, unlike the multicomponent approach, some attempts have been made in coupling the drug release and the hydrogel mechanics [127], treating the drug as a diffusant species, according to Table 3.

In 2008, an interesting work aiming to describe the transient analysis of temperature-sensitive neutral hydrogels was published [126] and improved in the years [128] until the last publication in 2012 [125]. The starting point was the mass and the momentum conservation of the phases, coupled with the energy conservation. It was demonstrated [126] that the energy transport is considerably quicker than the deformation; therefore, the system can be accurately described through an isothermal deformation based on the new equilibrium temperature [125].

Considering the conservation of the mass for the two phases (Table 3 equations (B)):

$\begin{array}{c}\frac{\partial {\varphi }_1}{\partial t}=-\nabla \cdot \left({\varphi }_1{v}_1\right)\\ \frac{\partial {\varphi }_2}{\partial t}=-\nabla \cdot \left({\varphi }_2{v}_2\right)\end{array}$E71

where ${\varphi }_2$ and ${\varphi }_1$ are the polymer and the fluid volume fractions and $v_2$ and $v_1$ the polymer and fluid velocities. Summing these equations, the mixture mass conservation expression was derived:

$0=-\nabla \cdot \left({\varphi }_2{v}_2+\left(1-{\varphi }_2\right)v_1\right)=-\nabla \cdot v_{\mathrm{vol}}$E72

where $v_{\mathrm{vol}}$ is the system volume average velocity. After some manipulations, it is possible to reformulate the polymer mass conservation in terms of the relative velocity ($v_R=v_1-v_2$):

$\nabla \cdot v_2=-\frac{\nabla \cdot q}{{\rho }_{1,0}}$E73

where $q={\varphi }_1{\rho }_{1,0}{v}_R$ is the relative flow vector of the fluid phase with respect to the polymer phase.

Considering the momentum conservation for the polymer and fluid phases (Table 3 equations (H) for nonionic system and with $({\mu }_1\dot{\gamma })\to 0$):

$\begin{array}{c}0=\nabla \cdot \left(-{\varphi }_{1}p\delta \right)+f_{12}\\ 0=\nabla \cdot \left(-{\varphi }_{2}p\delta +{\sigma }_2\right)+f_{21}\end{array}$E74

In eq. 74, $\text{p}$ is the fluid pressure inside the hydrogel, ${\sigma }_2$ is the polymer stress tensor, and $f_{12}$ and $f_{21}$ are the reaction couple of the drag force between the two phases: $f_{12}=-f_{21}=p\nabla {\varphi }_P+\zeta v_R$, where $\zeta$ is the friction coefficient defined as $\zeta ={\mu }_1{\varphi }_1^2/ \kappa$, function of the fluid viscosity ${\mu }_1$, of the fluid volume fraction, and of the hydrogel permeability $\kappa$. From the momentum balance on the fluid phase, the pressure gradient can be related to $q$ :

$\nabla p=-\frac{{\mu }_1}{\kappa }\frac{q}{{\rho }_{1,0}}$E75

While summing both eq. 71, the mixture moment conservation can be obtained:

$0=\nabla \cdot \left(-pI+{\sigma }_2\right)$E76

Equation 76 highlights that in the multiphasic approach, due to the preliminary assumptions, the stress tensor is the one of the polymer network, while the solvent–polymer interaction $-$ and eventually the ion effect $-$ influences the momentum of the hydrogel through the pressure term.

3.3.5.1. Constitutive equations

The Helmholtz free energy of a hydrogel is constituted by the elasticity-related term and by the polymer–solvent mixing-related term. As usual, the elastic contribution is attributed to the polymer chain through the Cauchy stress, and the mixing term is considered through the osmotic pressure.

The mixing term of the Helmholtz free energy was considered through the presence of an osmotic pressure with the expression [129]

$p^{osm}=-\frac{k_{b}T}{V_m}\left[{\varphi }_P+\chi {\varphi }_P^2+\ln \left(1-{\varphi }_P\right)\right]$E77

where $\chi$ is the polymer–solvent interaction parameter, $k_b$ is the Boltzmann constant, $T$ is the temperature, and $V_m$ is the volume occupied by one monomer.

The Cauchy stress tensor for the polymer ${\sigma }_2$ was derived using eq. 46 and an affine elastic free energy:

${\sigma }_2=G\left(B-\frac{1}{2}\delta \right)$E78

Equation 78 is similar to eq. 64, a part from the factor 0.5 that multiplies the unit tensor. It is worth highlighting that in some cases, it might be useful to separate the concept of reference configuration (state of polymerization of the hydrogel) from the initial configuration to avoid singular behavior (i.e., for polymerization in gas where ${\varphi }_2\approx 1$, the osmotic pressure goes to infinity). However, this is not the case, and the reader is referred to [105, 119, 129] for further information.

At this point with the proper initial and boundary conditions, knowing the constitutive equations, the partial differential eqs. 73, 75, and 76 could be simultaneously solved for $q, p$, and $v_2$ that in conjunction with a moving mesh and the other algebraic equations previously defined, allow to completely characterize the hydrogel deformation due to a temperature variation. The authors, to avoid the use of a moving mesh, recast the equations in Lagrangian configuration; however, since this last step is not indispensable for the modeling approach comprehension, it will not show here, and the interested readers are referred to [125, 126].

4. Nomenclature

$A$	Generic area	$\left[m^2\right]$
$A$	Helmholtz free energy per unit volume	$\left[J/m^3\right]$
$A_{el}$	Elastic contribution to the Helmholtz free energy	$\left[J/m^3\right]$
$A_{ion}$	Ionic contribution to the Helmholtz free energy	$\left[J/m^3\right]$
$A_{mix}$	Mixing contribution to the Helmholtz free energy	$\left[J/m^3\right]$
$B$	Finger tensor (left Cauchy–Green)	$[-]$
$C$	Solvent concentration in the reference configuration (see 3.3.4)	$\left[1/m^3\right]$
$C_n$	Flory’s characteristic ratio	$[-]$
$c^{\text{*}}$	Solvent concentration at the glass–rubber interface (see 2.4)	$\left[k{g}_{solvent}/k{g}_{dryPolymer}\right]$
$c_i$	Concentration of the i^th species	$\left[mol/m^3\right]$
$c_x$	Cross-linking concentration	$\left[mol/m^3\right]$
$D_i$	Diffusion coefficient of the i^th species	$\left[m^2/s\right]$
$D^{\left(i\right)}/Dt$	Substantial derivative following the i^th species velocity	$\left[1/s\right]$
$dF/ds$	Slope of the force–penetration curve	$\left[N/m\right]$
$F$	Faraday’s constant	$\left[C/mol\right]$
$F$	Fickian contribution to drug release	$[-]$
$F$	Deformation gradient tensor	$[-]$
$f$	Frictional draw coefficient	$\left[\left(J\cdot s\right)/m^2\right]$
$f_{\mathrm{ij}}$	Drag force per unit volume of the phase $i$ on the phase $j$	$\left[N/m^3\right]$
$G$	Hydrogel shear modulus	$\left[Pa\right]$
$G_0$	Hydrogel initial shear modulus at the stress-free state (gelation condition)	$\left[Pa\right]$
$H$	Tablet height	$\left[m\right]$
$H$	Inverse of the deformation gradient tensor	$[-]$
$I$	Light intensity	$[-]$
$I$	Ionic strength	$\left[mol/m^3\right]$
$I_0$	Light intensity of the dry tablet	$[-]$
$I_B$	First invariant of the Finger tensor	$[-]$
$II{I}_B$	Third invariant of the Finger tensor	$[-]$
$I_{max}$	Light intensity of the fully hydrated gel	$[-]$
$J$	Solvent flux in the reference configuration (see 3.3.4)	$\left[1/(m^{2}s\right)]$
$j_i$	Mass flux of the i^th species	$\left[kg/\left(m^{2}s\right)\right]$
$K_a$	Acid dissociation constant
$K_b$	Basic dissociation constant
$k$	Hydraulic permeability	$\left[m^2\right]$
$k$	Generic constant
$k_b$	Boltzmann’s constant	$\left[J/K\right]$
$k_{er}$	Erosion constant	$\left[kg/\left(m^{2}s\right)\right]$
$k_{eros}$	Erosion constant velocity	$\left[m/s\right]$
$k_{swe}$	Constant of swelling	$[-]$
$L_{Ch}$	Characteristic diffusion path length	$\left[m\right]$
$L_c$	Length of the polymeric chain	$\left[m\right]$
$l$	Length of the bond along the polymer backbone	$\left[m\right]$
$M$	Average molecular weight	$\left[kg/mol\right]$
$M_{\infty }$	Amount of drug released at infinite time	$\left[kg\right]$
$M_f$	Weight of the polymeric chain	$\left[kg\right]$
$M_i$	Molecular weight of the i^th species	$\left[kg/mol\right]$
$M_r$	Molecular weight of the repeating unit	$\left[kg/mol\right]$
$M_t$	Cumulative amount of drug released at time $t$	$\left[kg\right]$
${\overline{M}}_C$	Molecular weight between two consecutive cross-links	$\left[kg/mol\right]$
${\overline{M}}_n$	Number-average molecular weight of the polymer	$\left[kg/mol\right]$
$m_i$	Mass of the i^th species	$\left[kg\right]$
$N$	Number of links for chain	$[-]$
$N_A$	Avogadro constant	$\left[1/mol\right]$
$N_{De.D}$	Diffusional Deborah number	$[-]$
$N_{Sw.A}$	Swelling area number	$[-]$
$N_{Sw.I}$	Swelling interface number	$[-]$
$N_c$	Number of species	$[-]$
$n_i$	Mole of the i^th species	$\left[mol\right]$
$P$	First Piola–Kirchhoff stress tensor	$\left[Pa\right]$
$p$	Total intrinsic fluid pressure	$\left[Pa\right]$
$p\text{'}$	Total intrinsic fluid pressure minus osmotic pressure (eq. L)	$\left[Pa\right]$
$p^{osm}$	Osmotic pressure	$\left[Pa\right]$
$Q$	Volumetric swelling ratio	$[-]$
$q$	Relative flow vector of fluid phase with respect to the polymer phase	$\left[kg/\left(m^{2}s\right)\right]$
$R$	Gas constant	$\left[J/\left(Kmol\right)\right]$
$R$	Swelling contribution to drug release
$r_f$	Polymer fiber radius	$\left[m\right]$
$r_s$	Size of the diffusing solute	$\left[m\right]$
${\left({\overline{r}}_0^2\right)}^{1/2}$	Unperturbed root-mean-square of the end-to-end distance for polymer chains between two neighboring cross-links	$\left[m\right]$
$T$	Temperature	$\left[K\right]$
$\text{T}2$	Spin–spin relaxation time	$\left[s\right]$
$T_{exp}$	Experimental temperature	$\left[K\right]$
$T_g$	Glass transition temperature	$\left[K\right]$
$t_i$	Stress acting on the gel surface	$\left[Pa\right]$
$V$	System volume	$\left[m^3\right]$
$V_0$	Volume of the dry polymer	$\left[m^3\right]$
$V_1$	Molar volume of the solvent	$\left[m^3/mol\right]$
$V_G$	Volume of the swollen gel	$\left[m^3\right]$
$V_P$	Volume of the dry polymer	$\left[m^3\right]$
$v$	Velocity of the solvent penetration	$\left[m/s\right]$
$v_f$	Total free volume	$\left[m^3\right]$
$v_{f.p}$	Polymer free volume	$\left[m^3\right]$
$v_{f.w}$	Water free volume	$\left[m^3\right]$
$v_R$	Relative velocity of the fluid phase with respect to the polymer phase	$\left[m/s\right]$
$v_{\mathrm{eros}}$	Erosion velocity	$\left[m/s\right]$
$v_i$	Velocity of the i^th species	$\left[m/s\right]$
$v_{\mathrm{mix}}$	Mass average velocity of the mixture	$\left[m/s\right]$
$v_{\mathrm{swe}}$	Swelling velocity	$\left[m/s\right]$
$v_{\mathrm{vol}}$	System volume average velocity	$\left[m/s\right]$
$W$	Total work done to the hydrogel	$\left[J\right]$
$W_P$	Penetration work	$\left[J\right]$
$W^{aff}$	Strain energy density function from affine network theory	$\left[J/m^3\right]$
$\dot{W_V}$	Energy dissipation due to viscous deformation	$\left[W/m^3\right]$
$\dot{W_s}$	Energy dissipation due to solvent migration	$\left[W/m^3\right]$
$x_i$	Molar fraction of the i^th species	$[-]$
$z_i$	i^th species charge	$[-]$
Greek symbols
$\alpha$	Elongation ratio	$[-]$
$\alpha$	Obstruction theory parameter	$[-]$
${\alpha }_f$	Thermal linear expansion coefficient of the polymer	$\left[1/K\right]$
$\beta$	Contribution of the water to the expansion of the polymer	$\left[k{g}_{dryPolymer}/k{g}_{solvent}\right]$
${\beta }_i$	Diffusivity parameter for the i^th species	$[-]$
$\gamma$	Parameters used to relate the light intensity to the weight fraction	$[-]$
$\dot{\gamma }$	Rate of strain tensor	$\left[1/s\right]$
$\delta$	Thickness of the swollen region through which the solute diffuses (see 1.2.5)	$\left[m\right]$
$\delta$	infinitesimal thickness of the boundary element (see 3.3.2)	$$\left[m\right]$$ .
$\delta$	Identity tensor	$[-]$
$\epsilon$	Relative permittivity	$[-]$
${\epsilon }_0$	Vacuum permittivity	$\left[F/m\right]$
$\zeta$	Friction coefficient	$\left[\left(Pa\cdot s\right)/m^2\right]$
${\theta }_D$	Characteristic time for the diffusion of the solvent in the polymer	$\left[s\right]$
${\lambda }_m$	Characteristic stress-relaxation time	$\left[s\right]$
$\overline{\mu }$	Chemical potential of the solvent in the external medium	$\left[J\right]$
${\mu }_{el}$	Number of junctions in the network	$[-]$
${\mu }_1$	Fluid viscosity	$\left[Pa\cdot s\right]$
${\nu }_{el}$	Number of elastic chains	$[-]$
$\xi$	Gel mesh size	$\left[m\right]$
$\eta$	Viscosity	$\left[Pa\cdot s\right]$
$\pi$	Pi
${\rho }_i$	Density of the i^th species	$\left[kg/m^3\right]$
$\sigma$	Cauchy stress tensor	$\left[Pa\right]$
$\tau$	Tensile stress	$\left[Pa\right]$
$\overline{\upsilon }$	Polymer-specific volume	$\left[m^3/kg\right]$
$\text{Φ}$	Total Helmholtz free energy	$\left[J\right]$
${\varphi }_i$	Volume fraction of the i^thspecies	$[-]$
${\chi }_{12}$	Flory–Huggins solvent–polymer interaction parameter	$[-]$
$\psi$	Electric potential	$\left[V\right]$
$\mathrm{\Omega }$	Volume of each solvent molecule in the gel	$\left[m^3\right]$
${\omega }_{\text{i}}$	Mass fraction of the i^th species	$[-]$