The hydration dynamics of polyelectrolyte gels with applications to cell motility and drug delivery

Abstract

We combine the physics of gels with the hydrodynamics of two-phase fluids to construct a set of equations that describe the hydration dynamics of polyelectrolyte gels. We use the model to address three problems. First, we express the effective diffusion constants for neutral and charged spherically distributed gels in terms of microscopic parameters. Second, we use the model to describe the locomotion of nematode sperm cells. Finally, we describe the swelling dynamics of polyelectrolyte gels used for drug release.

Appendix

Appendix A: elasticity and functional derivatives

To define the stress tensor for the polymer mesh, Tanaka et al. (1973) originally assumed that the gel could be treated as an elastic solid with an elastic energy given by:

$$ F_{{\text{el}}} = \int {\left[ {\mu \left( {u_{ik} - \frac{1} {3}u_{jj} \delta _{ik} } \right)^2 + \frac{1} {2}Ku_{jj}^2 } \right]{\text{d}}^3 x} $$

(A1)

Here the components of the strain tensor are $u_{ik} = \left( {\partial u_i /\partial X_k + \partial u_k /\partial X_i } \right)/2$, μ is the shear modulus, and K is the bulk modulus. The first term in the integral is the energy associated with shear deformations of the solid and the second term is the energy associated with changes in volume. Moreover, they assumed that the solvent was stationary and set the viscous drag force acting on the gel polymer equal to the force obtained from the functional derivative of this energy (Tanaka et al. 1973). Here we introduce a more complete description of the free energy of the polyelectrolyte gel by taking into account not only the elastic energy but the mixing, counterion, and polymer interaction energies as well (English et al. 1996).

The energy we will use is dependent on two parameters, the metric determinant, $\sqrt g = \frac{{\phi _{{\text{init}}} }} {\phi }$, and the trace of the metric, ${\text{Tr}}\left( {g_{ij} } \right) = \frac{{\partial X_i }} {{\partial x_j }}\frac{{\partial X_i }} {{\partial x_j }}$. As such, the energetic cost for deformations of the gel can be written as:

$$ F = \iiint_V {\left( {\Theta \left( \phi \right) + \frac{\phi } {{\phi _{{\text{init}}} }}\Upsilon \left( {{\text{Tr}}\left( {g_{ij} } \right)} \right)} \right)}_{} {\text{d}}^3 X $$

(A2)

with Θ and $\Upsilon$ functions that describe the energetic cost for volume and volume/shear deformations, respectively, and the integral being taken over the volume of the gel, V. This appendix defines notation and shows how functional derivatives are calculated and used to derive the forces per volume that act on the polymer mesh.

We begin by defining the position of points in the polymer mesh through the vector X. The initial position of the polymer mesh at time t=0 is x _init=x(t=0). From this, we can define the displacement vector for these material points as u=X−x _init. The volume fraction of the polymer is obtained from this by:

$$ \frac{{\phi _{{{\text{init}}}} }} {\phi } = det{\left[ {\frac{{\partial {\user2{X}}}} {{\partial {\user2{x}}_{{{\text{init}}}} }}} \right]} = det{\left[ {\hat{I} + \frac{{\partial {\user2{u}}}} {{\partial {\user2{x}}_{{{\text{init}}}} }}} \right]} = {\sqrt g } $$

(A3)

with φ_init the initial polymer volume fraction and $\hat I$ the identity matrix. Using this notation, we now discuss how to derive the force from a given energy function. To better elucidate the method, we begin by considering only deformations that change the volume. The variation in the energy is taken with respect to the initial position vector, x _init, over the initial volume, V ₀, as the limits of this integration are not altered by the variation:

$$ \delta F = \iiint_{V_0 } {\frac{\partial } {{\partial \phi }}}\left[ {\frac{{\phi _{{\text{init}}} }} {\phi }\Theta \left( \phi \right)} \right]\delta \phi _{} {\text{d}}^3 x_{{\text{init}}} $$

(A4)

It can be shown (Onuki 1989) that:

$$ \delta \phi = - \phi \sum\limits_{l,m} {\left[ {\frac{{\partial x_{{\text{init}},m} }} {{\partial X_l }}} \right]} \left[ {\frac{\partial } {{\partial x_{{\text{init}},m} }}\delta X_l } \right] = - \phi \sum\limits_l {\left[ {\frac{\partial } {{\partial X_l }}\delta X_l } \right]} $$

(A5)

Substituting this expression into Eq. (A4) and transforming back to the displaced coordinates, we obtain:

$$ \delta F = \iiint_V {\left( { - \frac{{\phi ^2 }} {{\phi _{{\text{init}}} }}} \right)}\frac{\partial } {{\partial \phi }}\left[ {\frac{{\phi _{{\text{init}}} }} {\phi }\Theta \left( \phi \right)} \right]\sum\limits_l {\left[ {\frac{\partial } {{\partial X_l }}\delta X_l } \right]} \,{\text{d}}^3 X $$

(A6)

This is then integrated by parts to give:

$$ \delta F = - \iint_{\partial V} {\sum\limits_{l,m} {n_m } \Pi _{lm}^{\text{p}} \delta X_l }\,{\text{d}}A + \iiint_V {\sum\limits_{l,m} {\left[ {\frac{\partial } {{\partial X_m }}\Pi _{lm}^{\text{p}} } \right]} }\delta X_l \,{\text{d}}^3 X $$

(A7)

where n is the unit outward normal vector to the surface element dA, ∂V is the surface of the gel, and the isotropic polymer stress due to changes in volume is defined as:

$$ \Pi _{lm}^{\text{p}} = \left( {\phi \frac{{\partial \Theta }} {{\partial \phi }} - \Theta } \right)\delta _{lm} $$

(A8)

This also gives a force per volume:

$$ {\user2{f}}^{{\text{p}}} = - \nabla _{X} \cdot \Pi ^{{\text{p}}} $$

(A9)

where $\nabla _X$ is the gradient over the displaced coordinates, X.

In a similar fashion, we can find the stress from the trace term:

$$ \sigma _{{\text{tr}},ij} = - \frac{\phi } {{N_{\text{x}} V_m }}\sum\limits_l {\frac{{\partial X_i }} {{\partial x_l }}} \frac{{\partial X_j }} {{\partial x_l }} $$

(A10)

and we have used that:

$$ \Upsilon = \frac{{\phi _{{\text{init}}} }} {{N_{\text{x}} V_m }}{\text{Tr}}\left( {g_{ij} } \right) $$

(A11)

Combining these two terms, the total stress on the polymer is:

$$ \sigma _{ij}^{\text{p}} = \Pi _{ij}^{\text{p}} + \sigma _{{\text{sh}},ij} $$

(A12)

Appendix B: components of the gel free energy

The total free energy for the gel component can be broken up into six independent components arising from the elasticity of the polymer, the entropic mixing due to the fluid, interactions between polymer chains, the pressure due to the counter ions, Coulombic interactions, and a term that depends on gradients in the volume fraction as:

$$ F = F_{{\text{el}}} + F_{{\text{mix}}} + F_{int} + F_{{\text{ion}}} + F_{{\text{coulomb}}} + F_{{\text{grad}}} $$

(B1)

We choose an elastic free energy based on rubber elasticity (Flory 1953):

$$ F_{{\text{el}}} = \frac{{k_{\text{B}} T}} {{2N_{\text{x}} V_{\text{m}} }}\iiint_V {\phi \left[ {{\text{Tr}}\left( {g_{ij} } \right) - 3 - ln\left( {\frac{{\phi _0 }} {\phi }} \right)} \right]\,{\text{d}}^3 X} $$

(B2)

where φ is the polymer volume fraction, g _ij is the metric between the initial coordinates, x, and the final coordinates, X, k _B is Boltzmann’s constant, T the absolute temperature, N _x the number of monomers between crosslinks, φ₀ is a material parameter that sets the unstressed volume fraction, V _m is the equivalent volume occupied by one monomer, and V is the volume of the gel. From these definitions, the number of chains per volume is φ/Ν _x V _m. For the isotropic case:

$$ g_{ij} = \left( {\frac{{\phi _{{\text{init}}} }} {\phi }} \right)^{{2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-\nulldelimiterspace} 3}} \hat I $$

(B3)

and:

$$ F_{{\text{el}}} = \frac{{k_{\text{B}} T}} {{2N_{\text{x}} V_{\text{m}} }}\iiint_V {\phi \left[ {\left( {\frac{{\phi _0 }} {\phi }} \right) - 1 - \frac{1} {3}ln\left( {\frac{{\phi _0 }} {\phi }} \right)} \right]{\text{d}}^3 X} $$

(B4)

We choose the entropic mixing free energy defined by Flory theory (Flory 1953):

$$ F_{{\text{mix}}} = \frac{{k_{\text{B}} T}} {{V_{\text{m}} }}\iiint_V {\left( {1 - \phi } \right)ln\left( {1 - \phi } \right){\text{d}}^3 X} $$

(B5)

and the solvent–polymer interaction energy representing two-body interactions is given by:

$$ F_{{\text{int}}} = \frac{{k_{\text{B}} T}} {{V_{\text{m}} }}\iiint_V {\chi \phi \left( {1 - \phi } \right){\text{d}}^3 X} $$

(B6)

where χ is the Flory parameter that measures the strength of interaction between the polymer chains (Flory 1953). It is common to assume that χ is independent of volume. However, in principle, it can be a function of φ; here we assume a dependence:

$$ \chi = \chi _0 + \chi _1 \phi ^2 $$

(B7)

to approximate the glass–rubber transition for swelling gels.

The contribution to the energy from the counterions is obtained by assuming an ideal solution of charges with each species of charge labeled by the subscript i (English et al. 1996):

$$ F_{{\text{ion}}} = - k_{\text{B}} TN_{\text{A}} \iiint_V {\sum\limits_i {C_i ln\left( {C_i } \right)} (1 - \phi )\,{\text{d}}^3 X} $$

(B8)

where C _i is the molar concentration of the ith mobile counterion and N _A is Avogadro’s number. The factor of (1−φ) in the integral is because the ions only reside in the fluid component of the mixture. Likewise, C _i is the concentration with respect to the fluid volume, dV ^f, not the total volume, dV ^T.

The contribution to the free energy from Coulombic interactions is:

$$ F_{{\text{coulomb}}} = \iiint_V {\sum\limits_i {z_i e\Phi _{{\text{self}}} C_i } \,{\text{d}}^3 X} + \iiint_V {z_{\text{F}} e\Phi _{{\text{self}}} C_{\text{F}} \,{\text{d}}^3 X} $$

(B9)

where z _i is the valence of the ith mobile counterion, e is the electronic charge, φ_self is the self-consistent potential created by the polymer charges and the mobile ions, and C _F is the fixed charge density of the polymer.

For cases where large differences in the volume fraction might exist over small lengths, such as volume-phase transitions in gels (Sekimoto et al. 1989; Tanaka 1997), it is common to include a term that defines the energy required to maintain the sharp gradient. As in Landau–Ginzburg theory and Cahn–Hilliard theory (Cahn and Hilliard 1958), we assume the form of the energy to be:

$$ F_{{\text{grad}}} = \iiint_V {\frac{\gamma } {2}\nabla \phi \nabla \phi \,{\text{d}}^3 X} $$

(B10)

where γ is a material parameter that has units of energy/length. This type of energy has already been suggested for gels (Tanaka 1997).

Appendix C: the counterion swelling pressure of the gel

In Eq. (B8), the energy of the counterions is that of an ideal solution of charges. We assume that the diffusion of the ions is much quicker than the motions of the solvent and polymer. Therefore, we use the steady-state concentrations for the ions inside the gel. Because the counterions are confined within the solvent, the concentration is taken with respect to the solvent volume, dV ^f=(1−φ)dV ^T, where dV ^T is a total unit volume of the gel. If we assume that there is an effective charge per monomer, α, on the polymer, the fixed charge density of the polymer is equal to αN/(dV ^f) where N=φ/V _m is the number of monomers per volume and V _m is the volume of a monomer. From this, the fixed charge density of the polymer is:

$$ C_{\text{F}} = \frac{{\alpha \phi }} {{(1 - \phi )V_{\text{m}} }} $$

(C1)

Assuming a univalent salt in the solvent, the concentrations of the positive and negative ions, C ₊ and C ₋, respectively, can be related to the fixed charge density via the Donnan equilibrium (Overbeek 1956). Equating the chemical potentials across the boundary:

$$ \begin{gathered} C_ + = \frac{1} {2}\left( {\left( {4C_{\text{b}}^2 + C_{\text{F}}^2 } \right)^{1/2} - C_{\text{F}} } \right) \hfill \\ C_ - = \frac{1} {2}\left( {\left( {4C_{\text{b}}^2 + C_{\text{F}}^2 } \right)^{1/2} + C_{\text{F}} } \right) \hfill \\ \end{gathered} $$

(C2)

where C _b is the bulk molar concentration of positive (or negative) ions outside the gel. Therefore, the total molar ion concentration inside the gel, C _ion=C ₊+C ₋, can be related to the fixed charge concentration:

$$ C_{{\text{ion}}} = 2C_{\text{b}} \sqrt {1 + \left( {\frac{{\alpha \phi }} {{2\left( {1 - \phi } \right)N_{\text{A}} V_{\text{m}} C_{\text{b}} }}} \right)^2 } $$

(C3)

Note that in the limit of the neutral gel, C _ion=2C _b (2C _b is the total ion bath concentration), while in the limit of a highly charged gel, the concentration is (1/V _m)(αφ/(1−φ)). As $\phi \to 1$, this concentration diverges. We expect that the effective fixed charge per monomer decreases with volume fraction due to the condensation of charges onto the polymer as the gel dehydrates. Demanding that in the absence of solvent all the charges are condensed onto the polymer sets $\alpha \left( {\phi = 1} \right) = 0$. To satisfy this constraint, we choose α=α₀(1−φ), with α₀ a constant.

Using Eqs. (C1) and (C2), we can find the Coulombic contribution to the energy:

$$ \begin{gathered} F_{{\text{coulomb}}} = \iiint_V {e\Phi _{{\text{self}}} (C_ + - C_ - )\,{\text{d}}^3 X} + \iiint_V {z_{\text{F}} e\Phi _{{\text{self}}} C_{\text{F}} \,{\text{d}}^3 X} \\ = \iiint_V { - e\Phi _{{\text{self}}} C_{\text{F}} \,{\text{d}}^3 X} + \iiint_V {z_{\text{F}} e\Phi _{{\text{self}}} C_{\text{F}} \,{\text{d}}^3 X} \\ = 0 \\ \end{gathered} $$

(C4)

Therefore, the electrostatic contribution to the total free energy is zero. This is only valid at low ionic strengths and when the diffusion of the ions is much quicker than the dynamics of the polymer.

The free energy for the gel can now be written completely in terms of the polymer volume fraction and a shear term. Using the free energy given in Appendix B and Eq. (A8) and (A10) to calculate the polymer stress, we find:

$$ \left( {\frac{{V_{\text{m}} }} {{k_{\text{B}} T}}} \right)\Pi _{ij}^{\text{p}} = h(\phi )\delta _{ij} - \frac{\phi } {{N_{\text{x}} V_{\text{m}} }}\sum\limits_l {\frac{{\partial X_i }} {{\partial x_l }}} \frac{{\partial X_j }} {{\partial x_l }} $$

(C5)

where:

$$ h(\phi ) = - ln(1 - \phi ) - \phi - \chi \phi ^2 + N_{\text{A}} V_{\text{m}} \left( {C_{{\text{ion}}} - 2C_{\text{b}} } \right) + \frac{\phi } {{2N_{\text{x}} }} + \gamma \left( {\phi \nabla ^2 \phi - \frac{1} {2}\left( {\nabla \phi } \right)^2 } \right) $$

(C6)

Experiments have shown that, at low volume fractions, the gel shear modulus is much less than the bulk modulus (Tanaka and Fillmore 1979). If we neglect the shear component [which is the same as using the isotropic form of the elastic free energy (Eq. B4)] and the term that comes from the gradient squared piece in the energy, the polymer stress in Eq. (C5) reduces to:

$$ \left( {\frac{{V_{\text{m}} }} {{k_{\text{B}} T}}} \right)\Pi _{ij}^{\text{p}} = h(\phi )\delta _{ij} $$

(C7)

where:

$$ h(\phi ) = - ln(1 - \phi ) - \phi - \chi \phi ^2 + N_{\text{A}} V_{\text{m}} \left( {C_{{\text{ion}}} - 2C_{\text{b}} } \right) + \frac{{\phi _0 }} {{N_{\text{x}} }}\left( {\frac{1} {2}\left( {\frac{\phi } {{\phi _0 }}} \right) - \left( {\frac{\phi } {{\phi _0 }}} \right)^{\frac{1} {3}} } \right) $$

(C8)

is the dimensionless swelling pressure and we have assumed that the parameters χ, φ₀, α₀, N _x, and γ are constants. We have also subtracted 2N _A V _m C _b to account for the ion pressure outside the gel.

If we only assume a purely entropic energy as given by Eq. (B5), then:

$$ h(\phi ) = - ln(1 - \phi ) - \phi $$

(C9)

Plugging this into Eqs. (5) and (7), we find that:

$$ \zeta \left( {\frac{{\partial u}} {{\partial t}} - v} \right) = - \left( {\frac{{k_{\text{B}} T\phi }} {{V_{\text{m}} \left( {1 - \phi } \right)}}} \right)\nabla \phi $$

(C10)

For small φ, this should reproduce normal diffusion. Using that C _m≈φ/V _m is the concentration of gel monomer per fluid volume, we get:

$$ \zeta {\left( {\frac{{\partial {\user2{u}}}} {{\partial t}} - {\user2{v}}} \right)} = - k_{{\text{B}}} TC_{{\text{m}}} \nabla C_{{\text{m}}} $$

(C11)

For normal diffusion, the drag coefficient, ζ_D, is defined as the drag per particle. In our representation, ζ is the drag per total volume. The relation between these is $\zeta = \zeta _{\text{D}} C_{\text{m}}$, which shows that these equations reproduce pure diffusive behavior when only entropic effects are considered and we recover the Einstein relation for the diffusion coefficient, $D = {{k_{\text{B}} T} \mathord{\left/ {\vphantom {{k_{\text{B}} T} \zeta }} \right. \kern-\nulldelimiterspace} \zeta }_{\text{D}}$

Appendix D: linear approximation

The gel equations are a coupled, nonlinear set of equations. To gain insight into the dynamics that they describe, we linearize the isotropic expansion equations about an initial polymer volume fraction, φ_init, setting φ=φ_init(1−ε), and using the rescalings $u \to Lu$, with L a relevant length scale. We use Eqs. (3), (5), and (7), and expand to first order in ε to obtain:

$$ \begin{array}{*{20}l} {\varepsilon \hfill} & { \approx \hfill} & {{ - \nabla {\user2{u}}} \hfill} \\ {{{\user2{u}}_{t} } \hfill} & { = \hfill} & {{{\left( {A\nabla ^{2} {\user2{u}} + B\nabla \times {\left( {\nabla \times {\user2{u}}} \right)}} \right)}} \hfill} \\ \end{array} $$

(D1)

where:

$$ A = \left( {1 - \phi _{{\text{init}}} } \right)\left( {{{4\tilde \mu } \mathord{\left/ {\vphantom {{4\tilde \mu } {3 + \phi _{{\text{init}}} }}} \right. \kern-\nulldelimiterspace} {3 + \phi _{{\text{init}}} }}\left( {{{\partial h} \mathord{\left/ {\vphantom {{\partial h} {\partial \phi _{{\text{init}}} }}} \right. \kern-\nulldelimiterspace} {\partial \phi _{{\text{init}}} }}} \right)} \right) $$

(D2)

and:

$$ B = \left( {1 - \phi _{{\text{init}}} } \right)\left( {{{\tilde \mu } \mathord{\left/ {\vphantom {{\tilde \mu } {3 + \phi _{{\text{init}}} }}} \right. \kern-\nulldelimiterspace} {3 + \phi _{{\text{init}}} }}\left( {{{\partial h} \mathord{\left/ {\vphantom {{\partial h} {\partial \phi _{{\text{init}}} }}} \right. \kern-\nulldelimiterspace} {\partial \phi _{{\text{init}}} }}} \right)} \right) $$

(D3)

with:

$$ \tilde \mu = \left( {{{V_{\text{m}} } \mathord{\left/ {\vphantom {{V_{\text{m}} } {k_{\text{B}} T}}} \right. \kern-\nulldelimiterspace} {k_{\text{B}} T}}} \right)\mu $$

(D4)

h is the dimensionless swelling pressure that is derived in Appendix C, and ${{\partial h} \mathord{\left/ {\vphantom {{\partial h} {\partial \phi _{{\text{init}}} }}} \right. \kern-\nulldelimiterspace} {\partial \phi _{{\text{init}}} }}$ the derivative of h with respect to φ evaluated at φ=φ_init. The latter is equivalent to the equation derived by Tanaka et al. (1973) with bulk modulus:

$$ K = \left( {{{k_{\text{B}} T} \mathord{\left/ {\vphantom {{k_{\text{B}} T} {V_{\text{m}} }}} \right. \kern-\nulldelimiterspace} {V_{\text{m}} }}} \right)\left( {1 - \phi _{{\text{init}}} } \right)\phi _{{\text{init}}} \left( {{{\partial h} \mathord{\left/ {\vphantom {{\partial h} {\partial \phi }}} \right. \kern-\nulldelimiterspace} {\partial \phi }}} \right) $$

(D5)

This equation can now be used to calculate the rate of change in polymer volume fraction and the fluid velocity for swelling of a gel near some initial configuration φ_init. This can then be compared with the work that has already been done on spherical swelling which has neglected fluid velocity.

We begin by assuming that all variables are only functions of r, the radial coordinate, and that all vectors are only in the radial direction, $u = u(r)\hat r$ and $v = v(r)\hat r$. This allows the conservation equation to be integrated and the solvent velocity is obtained in terms of the polymer velocity:

$$ v(r) = \frac{{c(t)}} {{\left( {1 - \phi } \right)r^2 }} - \frac{{\phi u_t }} {{\left( {1 - \phi } \right)}} $$

(D6)

with c(t) a time-dependent constant. For a finite solution at r=0, c(t)=0, giving (1−φ)v+φu _t =0, i.e. the center of mass of the gel is stationary. Thus v>0 for hydrating gels; the fluid flows into the gel at a velocity equal to the ratio of gel volume to fluid volume. Using this result, the first-order polymer velocity equation is:

$$ u_t \approx \left( {\tilde \mu + \phi _{{\text{init}}} \frac{{\partial h}} {{\partial \phi _{{\text{init}}} }}} \right)\nabla _r^2 u $$

(D7)

with $\nabla _r$ being the radial gradient. This gives a diffusive dynamics:

$$ u_t = - D\nabla _r^2 u $$

(D8)

with a diffusion coefficient:

$$ D = - \frac{{k_{\text{B}} T}} {{\zeta V_{\text{m}} }}\left( {1 - \phi _{{\text{init}}} } \right)\left( {\tilde \mu + \phi _{{\text{init}}} \frac{{\partial h}} {{\partial \phi _{{\text{init}}} }}} \right) $$

(D9)

[The minus sign in the diffusion constant is due to the fact that as the gel swells, the displacement increases. In Eq. (D9), D>0.] If we assume a neutral gel, α=0, and utilize reasonable values of φ₀=0.1, V _m=10⁻²² cm³, χ=0.6, N _x=60 (English et al. 1996), μ=200 dyn/cm², and ζ=2×10¹¹ dyn s/cm⁴ (Tanaka and Fillmore 1979), we find D=6.8×10⁻⁷ cm²/s with φ_init=0.05 and D=4.7×10⁻⁷ cm²/s with φ_init=0.025, which are comparable to the diffusion constants previously observed for gels with these polymer volume fractions and the ratio D ₅/D _2.5=1.45 is what is experimentally observed (Tanaka et al. 1973). Using the continuity equation (Eq. 2) and the linearized equations, we obtain the effective diffusion equation for the polymer volume fraction:

$$ \varepsilon _t = D\nabla _r^2 \varepsilon $$

(D10)

When the gel/solvent interface is sufficiently permeable, water can flow into the gel to satisfy any imposed boundary condition on the gel. Using a zero stress boundary condition sets the polymer volume fraction at the interface to the unstressed volume fraction and the gel will relax to this unstressed state on a characteristic time scale given by τ≈L ²/D. This time decreases as the charge per monomer, α, is increased. Therefore, a polyelectrolyte gel swells faster than a neutral gel, as is observed experimentally.

To calculate the swelling radius versus time for the linear approximation of the spherically symmetric gel, we use the result derived in Tanaka and Fillmore (1979):

$$ \Delta a(t) = \left( {\frac{6} {{\pi ^2 }}} \right)\Delta a_0 \sum {n^{ - 2} } exp\left( {\frac{{ - Dn^2 t}} {{a^2 }}} \right) $$

(D11)

where Δa is the change in radius of the gel, Δa ₀ is the total change in radius, D is the diffusion coefficient found in Eq. (D9) and the sum is over n.

Appendix E: equations for the swelling gel

In this Appendix, we derive a dynamic equation for the polymer volume fraction, φ, in a Lagrangian frame of reference with respect to the initial coordinates. The coupled dynamic equations for the gel can be reduced to one equation by neglecting shear. Symmetry conditions permit this for certain geometries (e.g., spheres and slabs). For other geometries, the magnitude of the shear term is a few orders of magnitude smaller than that of the volume deformation term for a wide range of physically relevant parameters (Tanaka and Fillmore 1979). Moreover, we will see that neglecting shear greatly simplifies the analysis, yet fits the experimental data closely.

For the purposes of our calculations, it is easier to work in a reference frame fixed on the initial position of the gel, since this frame does not change in time. In this frame, the continuity equation for the polymer volume fraction does not contain convection, and, therefore:

$$ \left. {\partial _{t} } \right|_{x} \phi = - \phi \nabla _{X} \cdot {\user2{u}}_{t} $$

(E1)

where $\left. {\partial _t } \right|_x$ represents the partial derivative with respect to time at fixed material point x _init. Derivatives with respect to the initial position are given by:

$$ \frac{\partial } {{\partial X_i }} = \left( {\frac{{\partial x_{{\text{init}},j} }} {{\partial X_i }}} \right)\frac{\partial } {{\partial x_{{\text{init}},j} }} $$

(E2)

and we can define the metric as:

$$ g_{ij} = \frac{{\partial X_k }} {{\partial x_{{\text{init}},i} }}\frac{{\partial X_k }} {{\partial x_{{\text{init}},j} }} $$

(E3)

Substituting Eq. (5) into Eq. (2) and adding and subtracting the divergence of the polymer velocity, yields:

$$ \nabla _{X} \cdot {\left( {\frac{1} {\zeta }{\left( {\phi - 1} \right)}{\user2{f}}^{{\text{p}}} + {\user2{u}}_{t} } \right)} = 0 $$

(E4)

Plugging this into Eq. (E1) and using Eqs. (E2) and (E3) leads to a dynamic equation for φ:

$$ \left. {\partial _t } \right|_x \phi = - \frac{{\phi ^2 }} {{\phi _{{\text{init}}} }}\frac{\partial } {{\partial x_{{\text{init}},i} }}\left( {\frac{{\phi _{{\text{init}}} }} {{\zeta \phi }}\left( {1 - \phi } \right)f_i^{\text{p}} } \right) $$

(E5)

where $f_i^{\text{p}}$ is the ith component of the polymer force and it has been used that the metric determinant:

$$ g^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-\nulldelimiterspace} 2}} = {{\phi _{{\text{init}}} } \mathord{\left/ {\vphantom {{\phi _{{\text{init}}} } \phi }} \right. \kern-\nulldelimiterspace} \phi } $$

(E6)

For a spherically symmetric gel:

$$ \left( {\frac{R} {r}} \right)^2 \frac{{\partial R}} {{\partial r}} = \frac{{\phi _{{\text{init}}} }} {\phi } $$

(E7)

where R is the radial coordinate and r is the radius of the initial state.

If shear is neglected, and we assume isotropy:

$$ g^{ij} = \left( {\frac{\phi } {{\phi _{{\text{init}}} }}} \right)^{{2 \mathord{\left/ {\vphantom {2 d}} \right. \kern-\nulldelimiterspace} d}} \delta ^{ij} $$

(E8)

where d is the dimension of space and $\delta ^{ij}$ is the Kroencker delta function, then the polymer force can be written completely in terms of φ and ${{\partial \phi } \mathord{\left/ {\vphantom {{\partial \phi } {\partial x_{{\text{init}}} }}} \right. \kern-\nulldelimiterspace} {\partial x_{{\text{init}}} }}$. This simplifies the dynamics by giving a single nonlinear PDE for φ rather than the coupled set of equations (Eqs. 2, 3, 4, 5), and greatly simplifies the problem.

For the drug delivery swelling problem, Eq. (E5) was solved numerically on a spherically symmetric volume with initial volume fraction φ₀=0.95 everywhere, except at the boundary point where we seed the low-energy swollen state with φ=0.1. We used boundary conditions (see Eq. F1) with the drag coefficient discussed in Appendix G. The microscopic energy parameters were found from fits to the experimental data (see Appendix G). Equation (E5) was solved on the fixed initial coordinates using a variable time step method (MATLAB routine ode15s) and a forward explicit numerical scheme.

Appendix F: boundary conditions for a swelling gel

Taking the functional derivative of the free energy produces both an integral over the volume and surface terms [see, for example, Eq. (A8)]. The integrand of the volume integral gives the force per volume that acts on the gel. The integrand of the surface area integral provides the natural boundary conditions on the stress. In circumstances where there are no external forces acting, these boundary conditions apply. For the swelling gel, the natural boundary conditions are:

$$ \begin{array}{*{20}l} {{\sigma ^{{\text{p}}} \cdot {\user2{\hat{n}}}} \hfill} & { = \hfill} & {0 \hfill} \\ {{\gamma \phi \nabla \phi \cdot {\user2{\hat{n}}}} \hfill} & { = \hfill} & {0 \hfill} \\ \end{array} $$

(F1)

Tokita and Tanaka (1991) found that $\zeta \propto \phi _{}^{1.5}$ when φ<1; this is expected from Flory scaling theory since the drag coefficient scales like 1/ξ², where ξ is the length correlation for the gel. Scaling theory suggests that $\xi \propto \phi _{}^{ - 3/4}$. We use this dependence for the drag coefficient. We assume a volume fraction dependence on γ, and find the best fits to the data when $\gamma \propto \phi ^{ - 1}$

Appendix G: parameter fits from swelling experiments

Our model depends on a number of microscopic parameters that describe the material and physical properties of gels. The data of Budtova and Navard (1998) allow us to fit a number of the microscopic parameters in this model and enable us to fit the dynamic behavior of the swelling of dry gels (see Fig. 2 and Fig. 3). In fig. 2 of Budtova and Navard (1998), the swelling curves for spherical poly(0.75 sodium acrylate–0.25 acrylic acid) gels of differing radii are plotted. The bath ion concentration for this plot is C _b=0. Based on the mole fraction of the crosslinking density to that of the monomer, an estimate of the number of monomers per crosslink for these gels is N _x=1600 (Budtova and Navard 1998). If we assume that the volume of a monomer is ~0.1 nm³ (English et al. 1996), then the only parameters left to fit are χ, α₀, ζ, and φ₀. ζ sets the overall time scale for the process, and is found by scaling the computer time (time step interval) to the experimental time, giving ζ(φ=1)≈10¹³ dyn s/cm⁴, which is comparable to experiment (Tokita and Tanaka 1991). χ₀, α₀, and φ₀ were fit to provide a swelling ratio comparable to that observed in the experiments (Budtova and Navard 1998). χ₁ was set by fitting the experimental swelling data to the numerical results. We find χ=0.6+2.3φ², α₀=1.0, and φ₀=0.5. These values are reasonable compared to what has been measured previously (English et al. 1996).

Appendix H: model equations for cell crawling

We model the dynamics of the 1-D cytoskeletal strip of gel with the following four equations:

$$ \frac{{\phi _{{\text{init}}} }} {\phi } = 1 + \frac{{\partial u}} {{\partial x}} $$

(H1)

$$ \frac{\partial } {{\partial x}}\left( {\phi \frac{{\partial u}} {{\partial t}} + \left( {1 - \phi } \right)v} \right) = 0 $$

(H2)

$$ \zeta \left( {v - \frac{{\partial u}} {{\partial t}}} \right) = - \frac{{\partial p}} {{\partial x}} $$

(H3)

$$ \zeta _{{\text{ex}}} \frac{{\partial u}} {{\partial t}} + \zeta \left( {\frac{{\partial u}} {{\partial t}} - v} \right) = \frac{{\partial \sigma ^{\text{p}} }} {{\partial x}} $$

(H4)

Here Eqs. (H1) and (H2) are 1-D conservation equations for the polymer (Eq. 3) and the gel (Eq. 2), respectively. Equation (H3) is the 1-D force balance for the liquid (Eq. 4). Finally, Eq. (H4) is the 1-D gel force balance (Eq. 5), where we have added the effective viscous drag between the cytoskeleton and the surface; ζ_ex is the corresponding viscous drag coefficient. The gel stress is:

$$ \sigma ^{\text{p}} (\phi ) = (k_{\text{B}} T/V_{{\text{site}}} )h(\phi ) $$

(H5)

where the stress due to changes in volume, h(φ) is given by Eq. (C8). The shear terms are absent in 1-D in both the polymer and liquid force balance equations. We solve the problem on the moving boundary domain r(t)<x<f(t), where r(t) and f(t) are the coordinates of the rear and front of the cell, respectively. The rates of movement of the cell edges are given by the gel velocities at the edges, plus the rate of polymerization at the front and depolymerization at the rear, respectively. Thus the boundary conditions are:

$$ \frac{{{\text{d}}r}} {{{\text{d}}t}} = V_{\text{d}} + \frac{{\partial u}} {{\partial t}}|_{x = r(t)} ,{\text{ }}\frac{{{\text{d}}f}} {{{\text{d}}t}} = V_{\text{p}} (r,f) + \frac{{\partial u}} {{\partial t}}|_{x = f(t)} ,{\text{ }}V_{\text{p}} (r,f) = \frac{{V_{\text{p}}^0 }} {{f - r - l_1 }} $$

(H6)

Here V _d is the constant depolymerization rate. The polymerization (protrusion) rate at the front is the decreasing function of the cell length, l(t)=f(t)−r(t); l ₁ is a model parameter.

We assume that the pH gradient regulates the external drag coefficient so that ζ_ex=3ζ(1+0.8 tanh((x−r−l ₂)/l ₃)), where l ₂ and l ₃ are model parameters. (The internal drag coefficient, ζ, is given in Table 2.) The external drag coefficient is small and almost constant under the cell body but constant and large under the lamellipod. We also assume that the number of dimers between the adjacent crosslinks in the gel increases significantly from 25 at the front to 75 at the rear. We describe this effect quantitatively with a function similar to drag coefficient spatial behavior: N _x=50−25 tanh((x−r−l ₂)/l ₃). Finally, we reason that the equilibrium volume fraction of the gel scales as $\phi \propto N_{\text{x}}^{1/2}$. Indeed, the equilibrium size of a long flexible filament is proportional to the square root of its contour length (Landau et al. 1980). Thus, the number of MSP dimers in unit volume scales like the number of dimers between the adjacent crosslinks divided by the equilibrium distance between these crosslinks, $N_{\text{x}} /N_{\text{x}}^{1/2}$. Therefore, we assume that:

$$ \phi _0 = 0.005N_{\text{x}}^{1/2} $$

(H7)

where the coefficient is set such that the equilibrium volume fraction at the front is approximately 0.035. as is observed experimentally. For simplicity, we keep other model parameters constant; they are given in Table 2. In simulations, we used the length scale L=10 μm, the time scale:

$$ {{\zeta V_{{\text{site}}} L^2 } \mathord{\left/ {\vphantom {{\zeta V_{{\text{site}}} L^2 } {k_{\text{B}} T}}} \right. \kern-\nulldelimiterspace} {k_{\text{B}} T}} \approx 0.01\;{\text{s}} $$

(H8)

and the force scale:

$$ {{k_{\text{B}} TL^2 } \mathord{\left/ {\vphantom {{k_{\text{B}} TL^2 } {V_{{\text{site}}} }}} \right. \kern-\nulldelimiterspace} {V_{{\text{site}}} }} \approx 100\;{\text{dyn}} $$

(H9)

We also used additional parameters: l ₁=5 μm, l ₂=12 μm, l ₃=(1/15) μm, V _d=0.2 μm/s, $V_{\text{p}}^{\text{0}} {\text{ = 0}}{\text{.4 }}\mu {\text{m/s}}$. We non-dimensionalized the model equations using these scales.

Then, we used Eqs. (H2) and (H3) to express the fluid velocity in terms of the hydrostatic pressure gradient:

$$ v = - \frac{\phi } {\varsigma }\frac{{\partial p}} {{\partial x}} $$

(H10)

(This equation is valid up to a constant, but this constant is almost zero due to the conditions at the leading edge, where v=0, and the pressure gradient is exponentially small.) From Eqs. (H3) and (H4) we find:

$$ \frac{{\partial p}} {{\partial x}}{\text{ = }}\frac{{\partial \sigma ^{\text{p}} }} {{\partial x}} - \varsigma _{{\text{ex}}} \frac{{\partial u}} {{\partial t}} $$

(H11)

Substituting into Eq. (H4), we obtain:

$$ \left( {\zeta _{{\text{ex}}} \left( {1 - \phi } \right) + \zeta } \right)\frac{{\partial u}} {{\partial t}} = \left( {1 - \phi } \right)\frac{{\partial \sigma }} {{\partial x}} $$

(H12)

where:

$$ \sigma = \sigma (\phi ),{\text{ }}\phi = \frac{{\phi _{{\text{init}}} }} {{1 + \frac{{\partial u}} {{\partial x}}}} $$

(H13)

We solve the equations on a uniform grid using a forward explicit integration method. At each time step, we update the gel displacement, u, using the gel stress from the previous step. Then, we use Eqs. (H10) and (H11) to find the pressure gradient and liquid velocity. Then, using Eq. (H13) we compute the new values of the gel stress using zero stress boundary conditions. Finally, we use the boundary conditions (Eq. H6) to move the cell boundaries. We start from the initial conditions of zero stress and velocities everywhere and compute until the asymptotically stable steady stress, density, and velocity distributions evolve.