Also at ]Living Systems Materialogy Research Group, International Research
Frontiers Initiative, Tokyo Institute of Technology, 4259, Nagatsuta-cho,
Midori-ku, Yokohama 226-8501, Japan

# Simple model for the kinetics of stimuli-responsive gels with porous
structures

Yuhei Yamada<EMAIL_ADDRESS>Living Systems Materialogy Research
Group, International Research Frontiers Initiative, Tokyo Institute of
Technology, 4259 Nagatsuta-cho, Midori-ku, Yokohama 226-8501, Japan  Shingo
Maeda  [  Department of Mechanical Engineering, Tokyo Institute of Technology,
2-12-1 Ookayama Meguro-ku Tokyo 152-8550, Japan

###### Abstract

Stimuli-responsive gels are the gels that vary the volume depending on
environmental conditions. It has been reported that the stimuli-responsive gel
with an inhomogeneous structure exhibits faster volume change than the gel
without it. It is understood as a difference in the transfer dynamics of the
solvent, though, there are few models for discussing the effect of
inhomogeneity explicitly. In this paper, we propose a simple model for the
kinetics of volume change by introducing inhomogeneity as the probability
distribution for a random variable that characterizes the structure. Under
this framework, we demonstrate that inhomogeneity actually increases the rate
of volume change. Furthermore, through analysis of a simplified example, we
show that some features that are often seen in experiments can be reproduced
by the model. This model provides a simple equation for the time series of
volume change that can be used to characterize the experimental data based on
a simple physical picture of the phenomenon.

##  I Introduction

Stimuli-responsive gels are the gels that vary their volume depending on
environmental condition such as temperature, pH, irradiation, etck. Ahn  _et
al._ (2008). In recent years, they have attracted considerable attention in
engineering and biomedical fields as smart materials that adapt to the
environmentHines  _et al._ (2017); Lavrador  _et al._ (2021). In many
applications of the stimuli-responsive gels, there is demand of increasing
ratio and rate of volume change. The basic mechanism of the volume change is
explained by the Flory-Huggins theory, which is a statistical model that gives
the osmotic pressure of gels by calculating the free energy of the gel state
while considering the entropy and enthalpy for mixing polymer and solventDoi
(2013). The volume of gel is determined to balance the osmotic pressure and
elastic stress of polymer network. In the case of the stimuli-responsive gel,
the environmental condition changes the interaction between the polymer and
solvent, and it changes the osmotic pressure in the gel. Then, the gel absorbs
or expels the solvent, which results in a change in volume.

In this mechanism, the rate of volume change is limited by either the transfer
rate of stimulus (like heat) or the rate of polymer network’s motion. In most
cases, the motion of the polymer network is slower than the transfer of
stimulus; thus, it becomes a rate-limiting step of volume changeGehrke (1993).
The motion of a polymer network is well described by the Tanaka-Fillmore
theory, in which the motion of a polymer network follows a diffusion
processTanaka  _et al._ (1973); Li and Tanaka (1990); Doi (2009). The rate of
diffusion is characterized by the diffusion coefficient $D$, which equals to
the ratio of the elastic modulus of the polymer network to the friction
coefficient between the network and solvent. Many efforts are devoted to
enhance $D$ through control of the synthesis method and chemical conditions,
though, recent studies reveal another approach is very effective in increasing
the rate of volume change; introducing inhomogeneity. Experimental studies
revealed that when a responsive gel is made with a microporous open-celled
structure, the solvent can be absorbed or expelled by the convection process,
which makes the swelling/shrinking kinetics much faster than that by the
diffusion process of polymer networkKabra  _et al._ (1992); Takata  _et al._
(2002); Kato and Gehrke (2004); Gotoh  _et al._ (2004); Maeda  _et al._
(2008); Murase  _et al._ (2009); Seiffert (2012). It is also reported that
inhomogeneous structure drastically changes static properties of gels like
equilibrium swelling ratio and elasticityF. Di Lorenzo (2015); Lopez  _et al._
(2019). From theoretical perspective, while great development of numerical
simulation method enables to explore precise dynamics in the internal of
responsive gelsNikolov  _et al._ (2018); Camerin  _et al._ (2018), theoretical
understanding based on a simple model that captures essential features seems
to be less progressed. To understand the role of inhomogeneity in the
kinetics, we need a model that explicitly treats inhomogeneity. In this study,
we propose a simple model for the kinetics of responsive gels that introduces
inhomogeneity as a probability distribution. Under the framework, we first
discuss the role of the inhomogeneous structure and show that the
inhomogeneity actually makes the kinetics faster. Furthermore, we show that
some features that are often seen in experiments can be reproduced by the
model.

##  II Model

We focus on the uniform deformation of gels. Then, the volume change in
infinitesimal time variation is equal to the flow rate through the surface.
Denoting the volume of gel as $V$, the flow rate of solvent into the gel as
$Q$ and, time as $t$,

$\frac{dV}{dt}=Q.$ (1)

Our strategy is to model the flow of solvent in porous gels by the Hagen-
Poiseuille flow and introduce an inhomogeneous structure as the distribution
of the flow channel. In the absorbing and expelling process of the solvent in
actual gels, the solvent flows through tortuous channels. Here, we approximate
the channels as a bundle of capillaries with different radii. Denoting the
radius of capillary $j$ as $\lambda_{j}$ and flow rate of that as
$q(\lambda_{j})$, which is positive for the direction from external to
internal of the gel, the total flow rate into the gel is

$Q=\sum_{j}q(\lambda_{j}).$ (2)

For $q(\lambda_{j})$, we assume the Hagen-Poiseuille law, which is plausible
for flow in cylindrical pipes at low Reynolds numbers, as follows:

$q(\lambda_{j})=\gamma\frac{\lambda_{j}^{4}}{l}\Delta p,$ (3)

where $l$ is the length of the channel, $\Delta p$ is the pressure difference
between the external and internal of the gel, and $\gamma$ is a constant. The
approximation for $Q$ is similar to that used for modeling flow in solid
porous media, which is found to be valid for the casesWu and Yu (2007); Yu
(2008). $\Delta p$ is computed from various factors reflecting each mechanism
of stimuli-responsiveness and is given for individual cases. Although $l$ and
$\Delta p$ can be different for each channel, we assume them to be identical
for simplicity here. Note that $\gamma,\lambda_{j}$ and $l$ are positive from
their definition, though, $\Delta p$ can be positive and negative. The sing of
$\Delta p$ determines the direction of flow and so we can apply this framework
to both swelling and shrinking process.

For further simplification, we assume that the shape of the flow channel is
scaled by the total volume of channels (which equals to the volume of solvent)
as

$\lambda_{j}=A_{j}V_{s}^{\frac{1}{3}},\quad l=BV_{s}^{\frac{1}{3}},$ (4)

where $V_{s}$ is the volume of solvent and $B$ and $A_{j}$ are positive
constants. $A_{j}$ are random numbers and we denote the probability density
function as $p(A)$. $p(A)$ represents the inhomogeneity of the structure.
Denoting the volume of the polymer as $V_{p}$,

$\displaystyle V$ $\displaystyle=V_{p}+V_{s}$
$\displaystyle=V_{p}+\sum_{j}\pi\lambda_{j}^{2}l$ $\displaystyle=V_{p}+\pi
BV\sum_{j}A_{j}^{2}.$ (5)

When the number of $j$ is sufficiently large, the summation is replaced by an
integral. Then, from Eq. (5) we obtain

$\int A^{2}f(A)dA=\frac{1}{\pi B}\left(1-\frac{V_{p}}{V}\right).$ (6)

This becomes the normalization condition for $p(A)$ to assure that $V,V_{p}$
and $B$, i.e., the whole volume, volume of polymer and length of polymer, are
fixed. (Note that when $V$ and $V_{p}$ are given, $V_{s}=V-V_{p}$ is fixed and
so fixing $B$ is equivalent to fixing $l$.) Substituting Eqs. (2), (3), (4)
into Eq. (1) and replacing te summation with an integral, we obtain

$\displaystyle\frac{dV}{dt}$ $\displaystyle=\frac{\gamma}{B}V_{s}\int\Delta
pA^{4}f(A)dA$ $\displaystyle=\frac{\gamma}{B}(V-V_{p})\int\Delta
pA^{4}f(A)dA.$ (7)

When the volume fraction of polymer in the gel is sufficiently small, (which
is often valid for responsive gels,) $V_{p}$ becomes negligible compared with
$V$. Then, we obtain

$\frac{dV}{dt}=\frac{\gamma}{B}V\int\Delta pA^{4}f(A)dA.$ (8)

This becomes the differential equation for $V$ when $\Delta p$ and $f(A)$ are
given.

Now, we consider a case that $\Delta p$ is composed of osmotic pressure by
chemical concentration and resilience force by elasticity of structure. For
simplicity, We assume $\Delta p$ only depends on $V$ as

$\Delta p(V)=\frac{\alpha}{V}-\beta\left(\frac{V-V_{n}}{V_{n}}\right),$ (9)

where $\alpha,\beta$ and $V_{n}$ are constants that represent the strength of
osmotic pressure, strength of resilience force, and natural volume of the gel,
respectively. To consider physically meaningful cases, we assume $\beta>0$ and
$V_{n}>0$. Substituting Eq. (9) into Eq. (8), we obtain

$\frac{dV}{dt}=\gamma^{\prime}\left\\{\alpha-\beta\frac{V}{V_{n}}\left(V-V_{n}\right)\right\\},$
(10)

where

$\gamma^{\prime}\equiv\frac{\gamma}{B}\int A^{4}f(A)dA.$ (11)

This equation is equivalent to the logistic equation. Although there are four
parameters in Eq. (10), since $\alpha,\beta$ and $\gamma^{\prime}$ attribute
to the equation only in the form of $\alpha$ and $\gamma^{\prime}\beta$, the
number of relevant parameters is three. Figure 1 shows the schematics of the
phase space of Eq. (10). $V_{e}$ indicates the equilibrium volume
corresponding to a fixed point of the dynamics. Requiring the existence of a
stable fixed point, $\alpha>-\beta V_{n}/4$ must be assumed. When $\alpha<0$,
there is a solution that $V\to-\infty$. This is physically nonsense solution
due to the simplicity of the assumption for $\Delta p$. Thus, we ignore such
an initial condition in the following discussion.

Figure 1:  Phase space of the differential equation (10). Red arrows show the
direction of vector field.

##  III Analysis

With the initial condition ($V(0)=V_{0}$), we can integrate the differential
equation (10) as

$V(t)=\frac{-a}{ce^{t/\tau}+1}+V_{e},$ (12)

where

$\displaystyle\begin{split}a=V_{n}\left(1+\frac{4\alpha}{\beta
V_{n}}\right)^{\frac{1}{2}},\quad\tau=\frac{1}{\gamma^{\prime}\beta}\left(1+\frac{4\alpha}{\beta
V_{n}}\right)^{-\frac{1}{2}}\\\
V_{e}=V_{n}\left\\{\frac{1}{2}+\frac{1}{2}\left(1+\frac{4\alpha}{\beta
V_{n}}\right)^{\frac{1}{2}}\right\\},\quad
c=\frac{V_{e}+V_{0}-V_{n}}{V_{e}-V_{0}}.\end{split}$ (13)

The four parameters described by Eqs. (13) correspond to the set of the three
relevant parameters in Eq. (10) and the initial value of $V$. $V_{e}$ suffices
$\lim_{t\rightarrow\infty}V(t)=V_{e}$ as well as $\Delta p(V_{e})=0$. When
$\alpha/\beta V_{n}$ is small, $V_{e}\simeq V_{n}+\alpha/\beta$, which means
that the equilibrium volume is the sum of the natural volume and the ratio of
the strength of the osmotic pressure and the resilience force. $\tau$
characterizes the relaxation time of the exponential dependency, which
decreases when $\gamma^{\prime}$ and $\alpha$ becomes larger.

Figure 2 shows three typical behaviors of the model. Parameters are set as
$\alpha=0.1,\beta=10,\gamma^{\prime}=1,V_{n}=1$, and the case for three
initial conditions, $V_{0}=0.1,0.6,1.6$, are depicted. From Eq. (13),
$V_{e}\approx 1.1$ and $\tau\approx 0.1$ for these cases. Each curve
represents the typical behavior of large swelling, small swelling, and
shrinking. When $V_{0}$ is sufficiently small, $V(t)$ has an inflection point
at which the sign of curvature changes. When $V_{0}$ is not so small, $V(t)$
does not have an inflection point. From Fig. 1, it is found that the
inflection point is $V=V_{n}/2$, so the condition to have an inflection point
is $V_{0}/V_{n}<1/2$. When $c\gg 1$, the solution (12) becomes

$V(t)\simeq-\frac{a}{c}e^{-t/\tau}+V_{e},$ (14)

which indicates simple exponential relaxation.

Figure 2:  Typical time evolution of volume in the model.
($\alpha=0.1,\beta=10,\gamma^{\prime}=1,V_{n}=1$ and $V_{0}=0.1,0.6,1,6$.)

We now discuss the effect of inhomogeneity on the rate of volume change. In
our model, inhomogeneity is introduced as the randomness of $A_{j}$, which
characterizes the radius of cylindrical pores. Consider the case where the
probability density function is a Gaussian distribution:

$f(A)=\frac{C}{\sqrt{2\pi\sigma^{2}}}\exp{\left(-\frac{(A-\mu)^{2}}{2\sigma^{2}}\right)},$
(15)

where $\mu,\sigma$ and $C$ are constants that represent the mean, standard
deviation, and normalization constant, respectively. To ensure $A_{j}>0$, we
assume $\mu\gg\sigma$ here. Now, we compare the cases in which all physical
parameters, $V,V_{p},l,\alpha,\beta,\gamma,V_{n}$ and $\mu$, are the same
except $\sigma$. When $V,V_{p}$ and $l$ are fixed, we can compute $C$ by Eq.
(6) as

$C=\frac{1}{\pi B}\left(1-\frac{V_{p}}{V}\right)\frac{1}{\mu^{2}+\sigma^{2}}.$
(16)

Then, $\gamma^{\prime}$ is computed from Eq. (11) as

$\gamma^{\prime}=\frac{\gamma}{\pi
B^{2}}\left(1-\frac{V_{p}}{V}\right)\frac{1}{\mu^{2}+\sigma^{2}}(\mu^{4}+6\mu^{2}\sigma^{2}+3\sigma^{4}).$
(17)

It is found from Eq. (10) that the rate of volume change is proportional to
$\gamma^{\prime}$. Bacause the variation of $\sigma$ only affects
$\gamma^{\prime}$ in the right-hand side of Eq. (10), one can discuss the
effect of $\sigma$ through comparing the value of $\gamma^{\prime}$:

$\gamma^{\prime}(\sigma\neq 0)-\gamma^{\prime}(\sigma=0)=\frac{\gamma}{\pi
B^{2}}\left(1-\frac{V_{p}}{V}\right)\frac{5\mu^{2}\sigma^{2}+3\sigma^{4}}{\mu^{2}+\sigma^{2}}\geq
0.$ (18)

This indicates that the rate of volume change is larger for the inhomogeneous
case ($\sigma\neq 0$) than for the homogeneous case ($\sigma=0$).

##  IV Discussion

Although Eq. (12) is derived by assuming a very simple description for $\Delta
p$ as Eq. (9), it would be useful for characterizing the time series of volume
change obtained experimentally. It is often reported that the time series of
volume change, both swelling and shrinking process, is well fitted by
exponential function Li and Tanaka (1990); Suzuki and Hara (2001). Equation
(14), which is derived from Eq. (12) with $c\gg 1$, is consistent with the
behaviors. From Eqs. (13), $c\gg 1$ is equivalent to $V_{0}/V_{n}\gg 1$. This
suggests that the exponential relaxation of volume is observed when the
elastic resilience force sufficiently affects the kinetics. Although it should
be noted that the Tanaka-Fillmore theory can also reproduce such exponential
relaxation Li and Tanaka (1990).

When $V_{0}$ is sufficiently smaller than $V_{n}$, there emerges an inflection
point as discussed. This reflects the balance between the hydrodynamic and
elastic effect: in the swelling process, an increase in volume reduces the
hydrodynamic resistance of capillaries, which enhances the flow rate, whereas
approaching $V$ to $V_{e}$ suppresses the driving force, $\Delta p$. In the
model (Eq. (10)), the dominant effect changes at $V=V_{n}/2$. Such transition
of the dominant effect would be universal even when $\Delta p$ is a more
realistic complex function. Previous experimental studies of responsive gels
have often reported that the time series of volume change exhibits non-single
exponential behaviors such as a sigmoid curveKabra  _et al._ (1992); Kaneko
_et al._ (1995). These behaviors are explained through various ideas such as
formation of skin layer, domain structure in gel, etc., and they seem to be
reasonable for each case. However, here we suggest a possibility that such
behavior can be observed in some cases without considering the special
dynamics in the gel.

##  V Conclusion

In this paper, we propose a simple model for the volume change of stimuli-
responsive gels with a porous structure. We approximated the structure using a
bundle of capillaries with different radii and modeled the rate of volume
change to be equal to the hydrodynamic flow rate. In the framework, the
driving force is the pressure difference, $\Delta p$, between the internal and
external of the gel, which should be given for individual cases. For example,
we discussed the simplified case when $\Delta p$ consists of the osmotic
pressure and elastic resilience force. From the analysis, we found that the
inhomogeneity of the structure increases the rate of volume change.
Furthermore, we found that depending on the initial volume, the time series of
the volume exhibits the exponential relaxation or sigmoidal behavior. We
expect that the model will explain some of the experimental results on the
kinetics of stimuli-responsive gels.

## References

  * k. Ahn  _et al._ (2008) S. k. Ahn, R. M. Kasi, S.-C. Kim, N. Sharma, and Y. Zhou, Soft Matter 4, 1151 (2008).
  * Hines  _et al._ (2017) L. Hines, K. Petersen, G. Z. Lum, and M. Sitti, Advanced Materials 29, 1603483 (2017).
  * Lavrador  _et al._ (2021) P. Lavrador, M. R. Esteves, V. M. Gaspar, and J. F. Mano, Advanced Functional Materials 31, 2005941 (2021).
  * Doi (2013) M. Doi,  _Soft matter physics_ (Oxford University Press, 2013).
  * Gehrke (1993) S. H. Gehrke, Synthesis, equilibrium swelling, kinetics, permeability and applications of environmentally responsive gels, in  _Responsive Gels: Volume Transitions II_ , edited by K. Dušek (Springer Berlin Heidelberg, Berlin, Heidelberg, 1993) pp. 81–144.
  * Tanaka  _et al._ (1973) T. Tanaka, L. O. Hocker, and G. B. Benedek, The Journal of Chemical Physics 59, 5151 (1973).
  * Li and Tanaka (1990) Y. Li and T. Tanaka, The Journal of Chemical Physics 92, 1365 (1990).
  * Doi (2009) M. Doi, Journal of the Physical Society of Japan 78, 052001 (2009).
  * Kabra  _et al._ (1992) B. Kabra, M. Akhtar, and S. Gehrke, Polymer 33, 990 (1992).
  * Takata  _et al._ (2002) S. Takata, K. Suzuki, T. Norisuye, and M. Shibayama, Polymer 43, 3101 (2002).
  * Kato and Gehrke (2004) N. Kato and S. H. Gehrke, Colloids and Surfaces B: Biointerfaces 38, 191 (2004).
  * Gotoh  _et al._ (2004) T. Gotoh, Y. Maeda, Y. Nakatani, and S. Sakohara, Journal of Chemical Engineering of Japan 37, 597 (2004).
  * Maeda  _et al._ (2008) S. Maeda, Y. Hara, R. Yoshida, and S. Hashimoto, Angewandte Chemie International Edition 47, 6690 (2008).
  * Murase  _et al._ (2009) Y. Murase, S. Maeda, S. Hashimoto, and R. Yoshida, Langmuir 25, 483 (2009).
  * Seiffert (2012) S. Seiffert, Macromolecular Rapid Communications 33, 1135 (2012).
  * F. Di Lorenzo (2015) S. S. F. Di Lorenzo, Polym. Chem. 6, 5515 (2015).
  * Lopez  _et al._ (2019) C. G. Lopez, A. Scotti, M. Brugnoni, and W. Richtering, Macromolecular Chemistry and Physics 220, 1800421 (2019).
  * Nikolov  _et al._ (2018) S. Nikolov, A. Fernandez-Nieves, and A. Alexeev, Applied Mathematics and Mechanics 39, 47–62 (2018).
  * Camerin  _et al._ (2018) F. Camerin, N. Gnan, L. Rovigatti, and E. Zaccarelli, Scientific Reports 8, 14426 (2018).
  * Wu and Yu (2007) J. Wu and B. Yu, International Journal of Heat and Mass Transfer 50, 3925 (2007).
  * Yu (2008) B. Yu, Applied Mechanics Reviews 61, 050801 (2008).
  * Suzuki and Hara (2001) A. Suzuki and T. Hara, The Journal of Chemical Physics 114, 5012 (2001).
  * Kaneko  _et al._ (1995) Y. Kaneko, R. Yoshida, K. Sakai, Y. Sakurai, and T. Okano, Journal of Membrane Science 101, 13 (1995).