Stochastic model of randomly end-linked polymer network micro-regions

TL;DR

This paper develops a stochastic model for the formation of polymer networks, focusing on the variability of micro-region compositions during gelation, providing insights into network dynamics and heterogeneity.

Contribution

It introduces a master equation approach to model the stochastic formation and variability of micro-regions in polymer networks, including effects of annealing and cooperative binding.

Findings

Model captures the evolution of micro-region configurations over time.

Quantifies variability of network micro-regions at different reaction extents.

Highlights the impact of cooperative binding on network formation.

Abstract

Polymerization and formation of crosslinked polymer networks are important processes in manufacturing, materials fabrication, and in the case of hydrated polymer networks, synthesis of biomedical materials, drug delivery, and tissue engineering. While considerable research has been devoted to the modeling of polymer networks to determine averaged, mean-field, global properties, there are fewer studies that specifically examine the variance of the composition across "micro-regions" (composed of a large, but finite, number of polymer network strands) within the larger polymer network.Here, we mathematically model the stochastic formation of polymer networks comprised of linear homobifunctional network strands that undergo an end-linking gelation process. We introduce a master equation that describes the evolution of the probabilities of possible network micro-region configurations as a function of time and extent of reaction. We specifically focus on the dynamics of network formation and the statistical variability of the gel micro-regions, particularly at intermediate extents of reaction. We also consider possible annealing effects and study how cooperative binding between the two end-groups on a single network-strand affects network formation. Our results allow for a more detailed and thorough understanding of polymer network dynamics and variability of network properties.

Tables

Table 1. Table 1: Summary of variables used.

Symbol	Representation
$N$	Number of reactive end-groups per network strand
$\ell$	Number of bound end-groups per network strand
$s_{\ell }$	Designation of strand type with $\ell$ bound end-groups
$N_{\text{s}}$	Number network strands per micro-region
$n_0$	Number of $s_{\text{0}}\text{-network}$ strands per micro-region
$n_{\text{1}}$	Number of $s_{\text{1}}\text{-network}$ strands per micro-region
$n_2$	Number of $s_{\text{2}}\text{-network}$ strands per micro-region
$m$	Total number of bound end-groups per micro-region
$p$	Extent of reaction, $m/\left(N_{\text{s}}N\right)$
$t$	Time of reaction (time units)
$P(n_1,n_2,t)$	Micro-region configuration $\{n_0,n_1,n_2\}$ probability at $t$
$\alpha$	Reactivity/cooperative binding parameter (unitless)
$\lambda$	Binding rate (${\text{time}}^{-1}$)
$\kappa$	End-group rearrangement rate to binding rate ratio (unitless)

Equations

$$P_{\ell}(p)={N\choose\ell}{p}^{\ell}\big{(}1-p\big{)}^{N-\ell},$$

$$n_{\ell }(p)=N_{\text{s}}\left(\genfrac{}{}{0ex}{}{2}{\ell }\right)p^{\ell }(1-p{)}^{2-\ell }.$$

$$\mathcal{N}(n_0,n_1,n_2)=2^{n_1}\left(\genfrac{}{}{0ex}{}{N_{\text{s}}}{n_0{n}_1{n}_2}\right).$$

$$\mathcal{N}(N_{\text{s}},m,n_2)=\frac{2^{m-2n_2}{N}_{\text{s}}!}{(N_{\text{s}}-m+n_2)!(m-2n_2)!n_2!}.$$

$$Z_{N_{\text{s}},m}=\sum _{n_2=0}^{[m/2]}\mathcal{N}(N_{\text{s}},m,n_2),$$

$$P_{N_{\text{s}},m}(n_2)=\frac{\mathcal{N}(N_{\text{s}},m,n_2)}{Z_{N_{\text{s}},m}}.$$

$$⟨n_0⟩$$

$$⟨n_1⟩$$

$$⟨n_2⟩$$

$$\mathcal{N}(N_{\text{s}},m,n_2,\alpha )=\frac{(2/\alpha {)}^{m-2n_2}{N}_{\text{s}}!}{(N_{\text{s}}-m+n_2)!(m-2n_2)!n_2!},$$

$$P_{N_{\text{s}},m,\alpha }(n_2)=\frac{\mathcal{N}(N_{\text{s}},m,n_2,\alpha )}{\sum _{n_2=0}^{[m/2]}\mathcal{N}(N_{\text{s}},m,n_2,\alpha )}.$$

$$⟨n_2^2⟩=\sum _{n_2=0}^{[m/2]}{n}_2^2{P}_{N_{\text{s}},m,\alpha }(n_2),$$

$$\begin{array}{rl}\frac{\text{d}P(n_1,n_2,t)}{\text{d}t}=& 2\lambda (N_{\text{s}}-n_1-n_2+1)P(n_1-1,n_2,t)\\ & +\lambda \alpha (n_1+1)P(n_1+1,n_2-1,t)\\ & -\lambda [2(N_{\text{s}}-n_1-n_2)+\alpha n_1]P(n_1,n_2,t),\end{array}$$

$$\sum _{n_2=0}^{N_{\mathrm{s}}}\sum _{n_1=0}^{N_{\mathrm{s}}-n_2}1=\frac{(N_{\mathrm{s}}+2)(N_{\mathrm{s}}+1)}{2}.$$

$$\begin{array}{rl}\frac{\text{d}P(n_0,n_1,t)}{\text{d}t}=& 2\lambda (n_0+1)P(n_0+1,n_1-1,t)\\ & +\lambda \alpha (n_1+1)P(n_0,n_1+1,t)\\ & -\lambda \left(2n_0+\alpha n_1\right)P(n_0,n_1,t).\end{array}$$

$$⟨n_{\ell }(t)⟩=\sum _{n_1,n_2}{n}_{\ell }P(n_1,n_2,t),$$

$$\frac{\text{d}⟨n_0(t)⟩}{\text{d}t}=$$

$$\frac{\text{d}⟨n_1(t)⟩}{\text{d}t}=$$

$$\frac{\text{d}⟨n_2(t)⟩}{\text{d}t}=$$

$$⟨n_0(t)⟩=$$

$$⟨n_1(t)⟩=$$

$$⟨n_2(t)⟩=$$

$$⟨m(t)⟩=\frac{2N_{\text{s}}}{2-\alpha }\left(2-\alpha -e^{-\alpha t}+(\alpha -1)e^{-2t}\right).$$

$$⟨p(t)⟩=1-\frac{1}{2-\alpha }[e^{-\alpha t}-(\alpha -1)e^{-2t}],$$

$$⟨p(t)⟩=1-e^{-t}.$$

$$⟨n_0⟩=$$

$$⟨n_1⟩=$$

$$⟨n_2⟩=$$

$$G(z_0,z_1,t)=\sum _{n_0,n_1}P(n_0,n_1,t)z_0^{n_0}{z}_1^{n_1},$$

$$\frac{\partial G}{\partial t}=2(z_1-z_0)\frac{\partial G}{\partial z_0}+\alpha (1-z_1)\frac{\partial G}{\partial z_1}.$$

$$\begin{split}G(z_{0},z_{1},t)=&\Bigg{[}z_{0}e^{-2t}+\frac{2z_{1}\left(e^{-\alpha t}-e^{-2t}\right)}{2-\alpha}\\ &+1+\frac{\alpha e^{-2t}-2e^{-\alpha t}}{2-\alpha}\Bigg{]}^{N_{\text{s}}}.\end{split}$$

$$\begin{array}{rl}P(n_0,n_1,t)=& \left(\genfrac{}{}{0ex}{}{N_{\text{s}}}{n_0,n_1}\right)e^{-2t{n}_0}{\left(\frac{2(e^{-\alpha t}-e^{-2t})}{2-\alpha }\right)}^{n_1}\\ & \times {\left(1+\frac{\alpha e^{-2t}-2e^{-\alpha t}}{2-\alpha }\right)}^{N_{\mathrm{s}}-n_0-n_1}.\end{array}$$

$$\begin{array}{rl}P(n_1,n_2,t)=& \left(\genfrac{}{}{0ex}{}{N_{\text{s}}}{n_1,n_2}\right)e^{-2t(N_{\text{s}}-n_1-n_2)}\\ & \times {\left(\frac{2(e^{-\alpha t}-e^{-2t})}{2-\alpha }\right)}^{n_1}\\ & \times {\left(1+\frac{\alpha e^{-2t}-2e^{-\alpha t}}{2-\alpha }\right)}^{n_2}.\end{array}$$

$$⟨n_{\ell }^2(t)⟩=\sum _{n_1,n_2}{n}_{\ell }^{2}P(n_1,n_2,t)$$

Full text

I Introduction

The study of crosslinked polymer networks is important in many applications from heavy industry to biomedical research De Gennes (1979); Flory (1953); Rubinstein and Colby (2003); Peppas et al. (2006); Stauffer et al. (1982); Sangeetha and Maitra (2005); Laftah et al. (2011); Nicolson and Vogt (2001). Crosslinked polymer networks can be formed by various techniques, leading to a diverse and complex set of structures and properties. Of these network types, considerable attention has been paid to those formed by a process termed “end-linking”. End-linked networks are usually comprised of polymeric precursors, or “network strands,” that contain $N$ reactive end-groupsHild (1998); Mark and Sullivan (1977). During gelation, crosslinks, or “branchpoints,” link multiple end-groups together. For example, poly(ethylene glycol) (PEG)-based hydrogels, which are common in biomedical applications, are typically formed through the reaction of its end-groups Lin and Anseth (2009).

Network strands can bind to form a network via two main polymerization reaction mechanisms: step-growth and chain-growth. Several excellent reviews have been written to describe these processes Lin and Anseth (2009); STANFORD et al. (1984); Mespouille et al. (2009); Buwalda et al. (2014). Briefly, gelation by step-growth polymerization typically involves a defined binary reaction (e.g., thiol-ene or azide-alkyne reactions) between the network strand end-groups and the complementary binding sites of a multifunctional branchpoint which acts to crosslink the network strands together. Networks formed by step-growth polymerization are typically more homogeneous in structure since the number of functional groups per branchpoint can be fixed. Gelation by chain-growth polymerization occurs via a chain-extension reaction where the network strand end-groups bind to a growing chain of end-groups, termed the “active center” (e.g., free-radical polymerization of vinyl end-groups). The chain of end-groups forms a branchpoint that crosslinks the network strands together. Networks formed by chain-growth polymerization tend to have a more heterogeneous structure since the number of end-groups bound at the branchpoint is not fixed.

Network strands binding via either step- or chain-growth may exist in many states. For bifunctional strands with $N=2$ reaction sites three possibilities arise as depicted in Figure 1a: (i) the strand may be “free” where neither of the reactive ends have bound ($s_{0}$-strand); (ii) the strand may “dangle” where only a single end has bound and the strand dangles from the rest of the network ($s_{1}$-strand); or (iii) the strand may be “intact” where both ends are bound to the larger polymer network and bridge two different branchpoints ($s_{2}$-strand) Metters et al. (2000). Strands with both ends bound may also form a loop, where both ends are bound to the same branchpoint Tonelli and Helfand (1974). The proportion of free, dangling and intact network strands may affect the chemical and physical properties of the network, for example in water-swollen polymeric networks, bound strand ratios impact gel modulus, mesh size, and swellingRubinstein and Colby (2003).

Finally, the architecture of polymer networks formed by end-linking gelation is not spatially uniform. Heterogeneous domains within polymer networks exist that span a few to hundreds of nanometers in size and arise through variations in local strand concentration (termed “frozen concentration fluctuation”), heterogeneous distribution of crosslinking, or topological- and connectivity-based inhomogeneities due to variability in network strand assembly Gu et al. (2019); Seiffert (2017); Bastide and Leibler (1988); Ikkai and Shibayama (2005). We define these microscopic domains as “micro-regions.” For simplicity we assume micro-regions are statistically identical, independent, and composed of a fixed number $N_{\rm s}$ of strands. Each strand is also assumed to carry $N=2$ reactive end groups, the most representative experimental scenario Hild (1998), resulting in a total of $2N_{\text{s}}$ available binding sites per micro-region. We also denote by $m$ the number of bound end-groups in each micro-region, so that the fraction $p$ of bound end-groups per micro-region is $p=m/(2N_{\text{s}})$. This quantity is also known as the extent of reaction, and can be experimentally tuned to control the elastic modulus, viscosity, swelling, mesh size, and other network propertiesFlory (1953). By definition $0\leq p\leq 1$ since the number of bound end groups $m$ cannot exceed the total number of available ones $2N_{\text{s}}$. Finally we assume that micro-regions are large enough that boundary effects between bordering domains are negligible so that the free, dangling, intact strand distribution of two adjacent micro-regions are not correlated. Note that the same value of $0<p<1$, may be associated to different $\{n_{0},n_{1},n_{2}\}$ micro-region configurations with $n_{0}$ free, $n_{1}$ dangling, and $n_{2}$ intact strands. Figure 1b shows four different micro-region realizations within a larger network where $N_{\rm s}$ is fixed but different $\{n_{0},n_{1},n_{2}\}$ configurations arise, resulting in different extents of reaction $p$. In Figure 1d-f we show several distinct $\{n_{0},n_{1},n_{2}\}$ configurations corresponding to fixed $N_{\text{s}}$ and $p$.

Some quantities of interest may be derived using $p$ such as the likelihood P_{\ell}\big{(}p\big{)} De Gennes (1979); Flory (1953); Bibbo and Valles (1982); Dušek et al. (2002); Metters et al. (2000) of finding free ($\ell=0$), dangling $(\ell=1)$, and intact ($\ell=2$) strands for a given $p$. A stochastic analysis however would lead to an expression for the probability distribution of finding any micro-region configuration $\{n_{0},n_{1},n_{2}\}$ corresponding to a given $p$, offering a much richer understanding of the binding process. Previously developed stochastic models use a subset-of-states approachStanford and Stepto (1975); Stanford et al. (1975); Ahmad and Stepto (1980), where a polymerizing mixture is described as a set of “subgraph” states of monomeric strands, a subset of which is used to drive polymerizationSchamboeck et al. (2019); Wang et al. (2016). These models, however, only examine the connectivity of small subgraphs, typically made of only a few network strands, to represent large scale networks and predict bulk quantities such as the network gel point Lin et al. (2018); Schamboeck et al. (2019). Studies involving larger subgraphs containing a sizable number of strands (say, greater than ten) are still lacking. Finally, although several Monte Carlo numerical studies have examined network heterogeneityKroll and Croll (2017, 2015); Balabanyan et al. (2005); Gilra et al. (2000); Leung and Eichinger (1984a, b); Hosono et al. (2007), none of them have evaluated configuration probability distributions.

We aim to determine the probability distribution for a given configuration of free, dangling, and intact strands within a micro-region of $N_{\text{s}}$ bifunctional strands, given a total number of $m$ bound end-groups, or equivalently for fixed $p=m/2N_{\text{s}}$. This will allow us to go beyond the characterization of a micro-region by means of $p$ and $N_{\text{s}}$ alone and to obtain analytical expressions for micro-region properties that depend on possible $\{n_{0},n_{1},n_{2}\}$ configurations. In some experimental scenarios $m$ or $p$ may change among micro-region realizations, and it may be useful to understand the structure of the network for a given average, intermediate extent of reaction $\langle p\rangle$. Our results will thus be presented both as a function of time, and of the average extent of reaction. We draw on existing stochastic self assembly and nucleation modelsChou and D’Orsogna (2011); D’Orsogna et al. (2012); Yvinec et al. (2012); D’Orsogna et al. (2013); Davis and Sindi (2016) and utilize a master equation approach. Different forms of the master equation will be developed and analyzed to account for different end-group reactivities and the possibility for end-groups to dynamically rearrange within the micro-region. We do not model branchpoint functionality but focus on the number of intact, dangling, and free $\{n_{0},n_{1},n_{2}\}$ strands within micro-regions. As a result, the total number of branchpoints and topology of the network do not affect our modeling, so while the structures depicted in Figure 1 resemble those formed by chain-growth polymerization, our methods can be easily applied to step-growth polymerization as well. Table 1 lists the various quantities used in the remainder of this work.

II Mathematical Models and Analysis

For completeness, we first review basic combinatoric and equilibrium models describing network formation before introducing our master equation models.

II.1 Combinatoric and equilibium models

Many mathematical studies of network formation via end-linking have used combinatoric approaches to quantify the number of polymeric strands in a given state Miller and Macosko (1976); Macosko and Miller (1976); Metters et al. (2000); Tibbitt et al. (2013); Norris et al. (2017). The extent of reaction $p$, defined as the fraction of bound end-groups per micro-region, can also be interpreted as the probability that any end-group within a micro-region has bound. The probability $P_{\ell}(p)$ of finding an $s_{\ell}$-strand with $0\leq\ell\leq N$ bound end-groups is thus

[TABLE]

which assumes that of $N$ end-groups, $\ell$ are bound and $N-\ell$ are not. Equation 1 provides a basis for mean field end-linking gelation models used to predict network properties. Henceforth we assume $N=2$. The average number $n_{\ell}(p)$ of $s_{\ell}$-strands within a micro-region of $N_{\rm s}$ strands can thus be written as

[TABLE]

Equation 2 does not provide any information on the possible spatial arrangement of free, dangling, intact, strands within a micro-region. Some end-groups may also bind differently than others depending on their state. For example free strands might more readily bind than dangling ones since diffusion allows them to more easily navigate the local environment to find an appropriate reaction site. Cooperative binding arises when the unbound end-group of a dangling strand more readily binds to form a fully bound, intact strand due to its proximity to the polymerizing network, especially when the polymer solution is dilute. Uncooperative binding emerges when $s_{\text{2}}\text{-strand}$ formation from the binding of an existing $s_{\text{1}}\text{-strand}$ is hindered by negative allosteric effects, which has been shown to occur in rigid strands Hosono et al. (2007). Finally, the reaction steps associated with network formation may also be irreversible or reversible. Irreversible reactions lead to “quenched” network formation whose properties are highly dependent on initial conditions, while reversible reactions allow the network to rearrange while forming and “anneal.”

The binding scenarios described above lead to different probability distributions for a given micro-region configuration. Evaluating these distributions requires a more complex mathematical representation than Equations 1 or 2. Some can be can be derived via combinatoric arguments, for example in the case of reversible binding, when equilibrium is reached and the annealing process is complete. We evaluate such limit here and find the probability distribution for a given micro-region configuration $\{n_{0},n_{1},n_{2}\}$ with $n_{0}$ unbound $s_{\text{0}}\text{-strands}$, $n_{1}$ singly bound $s_{\text{1}}\text{-strands}$, and $n_{2}$ doubly bound $s_{\text{2}}\text{-strands}$, under the assumption that a total of $m$ end-groups have bound.

At equilibrium, the time and order at which strands were linked do not affect configuration likelihoods, so the task of finding the probability distribution for $\{n_{0},n_{1},n_{2}\}$ is equivalent to finding the number of ways ${\cal N}(n_{0},n_{1},n_{2})$ one can distribute $\{n_{0},n_{1},n_{2}\}$ among $N_{\text{s}}$ strands with $m$ total bindings. The above quantities are related by $n_{0}+n_{1}+n_{2}=N_{\text{s}}$ since all strands must be accounted for, and by $n_{1}+2n_{2}=m$ to include the contribution of each strand type to the total bound end-group count. Hence, a given micro-region with configuration $\{n_{0},n_{1},n_{2}\}$ can be equivalently described by $\{N_{\text{s}},m,n_{2}\}$. The extent of reaction $p$ can also be determined from $\{N_{\text{s}},m,n_{2}\}$ via $p=m/(2N_{\text{s}})$. Combinatoric arguments yield ${\cal N}(n_{0},n_{1},n_{2})$ as

[TABLE]

Here, the $2^{n_{1}}$ factor arises from the fact that the bound end-group on an $s_{\text{1}}\text{-strand}$ can be arranged in two configurations per strand. The above can be rewritten using $n_{0}=N_{\text{s}}-m+n_{2}$ and $n_{1}=m-2n_{2}$ as follows

[TABLE]

Upon summing over $n_{2}$ with $N_{\text{s}},m$ fixed, we find $Z_{N_{\text{s}},m}$ the partition function over all possible configurations, with $N_{\text{s}},m$ fixed

[TABLE]

where $[\cdot]$ indicates the integer part of its argument. The equilibrium probability distribution can finally be calculated as

[TABLE]

Equation 6 may be used to evaluate many different micro-region properties, such as averages, variances, and higher moments. We begin with the average number of free, dangling, and intact strands, respectively given by

[TABLE]

In Equations 7 the average, denoted by $\langle\cdot\rangle$, is taken across all micro-regions with the same $N_{s}$ and $m$, or equivalently, using all possible configurations within a single micro-region with $N_{\text{s}}$ strands and $m$ total number of bound end-groups. The above combinatoric argument assumes that end-group binding is accompanied by end-group annealing until equilibrium is reached, independent of the number of bound end-groups already present. However, within cooperative or uncooperative binding, bound end-groups may promote or hinder the binding of other end-groups. We include these phenomena by re-writing Equation 4 as

[TABLE]

where the reactivity parameter $\alpha>1$ represents cooperative binding, penalizing dangling ends in favor of $s_{\text{1}}\rightarrow s_{\text{2}}$ events. Values of $\alpha<1$ represent uncooperative binding where $s_{\text{0}}\rightarrow s_{\text{1}}$ events are favored. The neutral case is $\alpha=1$. Finally, the equilibrium probability distribution $P_{N_{\text{s}},m,\alpha}(n_{2})$ can be written as

[TABLE]

We plot $P_{N_{\text{s}},m,\alpha}(n_{2})$ in Equation 9 for several values of $n_{2}$, under three choices of $\alpha$ and as a function of the extent of reaction $p=m/(2N_{\text{s}})$ in Figures 2a–c. The solid lines connect micro-region configurations with the same $n_{2}$; we choose this representation as the number of intact “elastically effective” $s_{\text{2}}\text{-network}$ strands is an important feature of polymer networks and determines both the mechanical modulus and swelling behavior of the network Rubinstein and Colby (2003).

As $\alpha$ increases, all curves tend to shift to the left, as might be expected since increasing cooperative effects favor the emergence of $s_{2}$-strands for a given $p$. In Figure 2d we plot the average strand fractions $\langle n_{\ell}\rangle/N_{\text{s}}$ for $\ell=0,1,2$ as evaluated via Equations 7 for $N_{\text{s}}=40$ and as a function of $p$. Note that for any $\alpha$, the average quantity $\langle n_{1}\rangle$ is a symmetric function of $m$ about $N_{\rm s}$ as can be verified by imposing $m^{\prime}=2N_{\rm s}-m$ in Equation 7b and verifying that $\langle n_{1}\rangle$ remains unchanged. Since $p=m/2N_{\rm s}$, this also implies that $\langle n_{1}\rangle$ will be symmetric about $p=1/2$ for all values of $\alpha$, as seen in Figures 2e-f. We also calculate the second moment $\langle n_{2}^{2}\rangle$ defined as

[TABLE]

from which we obtain the variance $\mathrm{Var}(n_{2})=\langle n_{2}^{2}\rangle-\langle n_{2}\rangle^{2}$ where $\langle n_{2}\rangle$ is derived in Equation 7c. Similarly as for $\langle n_{1}\rangle$ one can verify that $\mathrm{Var}(n_{2})$ is symmetric about $p=1/2$ for all values of $\alpha$. Since $n_{1}=m-2n_{2}$, $n_{0}=N_{\text{s}}-m+n_{2}$, and given $\langle n_{1}\rangle$ and $\langle n_{0}\rangle$ from Equations 7a-b, $\mathrm{Var}(n_{1})=\langle n_{1}^{2}\rangle-\langle n_{1}\rangle^{2}$ and $\mathrm{Var}(n_{0})=\langle n_{0}^{2}\rangle-\langle n_{0}\rangle^{2}$ can also be derived using Equations 7c and 10. Figure 2e shows $\mathrm{Var}(n_{2})/N_{\rm s}^{2}$ as a function of $p$ for different values of $\alpha$. In each case, the maximum variance occurs when half of all possible end-groups have bound at $p=1/2$. As $\alpha$ deviates from the neutral condition $\alpha=1$, the bias towards certain bond types induced by cooperativity or uncooperativity causes the variance to decrease. In Figure 2f we plot $\mathrm{Var}(n_{2})/N^{2}_{\text{s}}$ as a function of $p$ for different values of $N_{\text{s}}$: the curve remains symmetric about $p=0.5$ and as $N_{s}$ increases, the normalized variance decreases.

II.2 Dynamic models: Master Equation approaches

We now derive the probability distribution $P(n_{0},n_{1},n_{2},t)$ of finding a given $\{n_{0},n_{1},n_{2},t\}$ micro-region configuration at time $t$ through a master equation that allows for the inclusion of reversible/irreversible (annealed/quenched) bond formation, and cooperative/uncooperative binding. Since the total number of strands per micro-region is constant, the constraint $n_{0}+n_{1}+n_{2}=N_{\text{s}}$ is obeyed at all times and effectively $P(n_{0},n_{1},n_{2},t)\rightarrow P(n_{1},n_{2},t)$. We compare equilibrium or steady state solutions to Equation 1; where possible we also determine the full time-dependent solution for $P(n_{1},n_{2},t)$, which can be used to derive other quantities of interest, such as the variance and higher moments.

II.2.1 Quenched end-group binding

The first case we consider is that of irreversible (or quenched) end-group binding, whereby once an end-group has bound, it will not detach. We also assume the binding rate $\lambda$ of an end-group is constant. Under these conditions, the master equation for the probability distribution $P(n_{1},n_{2},t)$ evolves according to

[TABLE]

where we have explicitly used the $n_{0}=N_{\text{s}}-n_{1}-n_{2}$ constraint. Equation 11 also includes the reactivity parameter $\alpha$: $\alpha>1$ represents cooperative binding so that $s_{\text{1}}\rightarrow s_{\text{2}}$ binding events are more likely than $s_{\text{0}}\rightarrow s_{\text{1}}$ events; the reverse is true for uncooperative binding, $\alpha<1$, where $s_{\text{0}}\rightarrow s_{\text{1}}$ events are favored. The first term on the right hand side of Equation 11 represents the process of an unbound strand attaching to the network structure to form a singly bound dangling strand ($s_{\text{0}}\rightarrow s_{\text{1}}$), which gives the configuration transition $\{n_{0}+1,n_{1}-1,n_{2}\}\to\{n_{0},n_{1},n_{2}\}$ (Figure 3).

The multiplicative factor $N_{\text{s}}-n_{1}-n_{2}+1$ represents the number $s_{\text{0}}\text{-strands}$ in the starting configuration that can bind to the network; the two prefactor is included since an $s_{\text{0}}\text{-strand}$ can bind to the network at either of its two unbound end-groups. Similarly, the second term represents an unbound end-group from a singly bound strand binding to the network and forming a doubly bound strand ($s_{\text{1}}\rightarrow s_{\text{2}}$). The related transition is $\{n_{0},n_{1}+1,n_{2}-1\}\to\{n_{0},n_{1},n_{2}\}$ (Figure 3). The multiplicative factor $n_{1}+1$ represents the number of $s_{\text{1}}\text{-strands}$ that can bind to the network to form an $s_{\text{2}}\text{-strand}$. Finally the last term describes the processes that drives the system out of the $\{n_{0},n_{1},n_{2}\}$ configuration, where either an $s_{\text{0}}\rightarrow s_{\text{1}}$ transition, with $\{n_{0},n_{1},n_{2}\}\to\{n_{0}-1,n_{1}+1,n_{2}\}$, or an $s_{\text{1}}\rightarrow s_{\text{2}}$ transition, with $\{n_{0},n_{1},n_{2}\}\to\{n_{0},n_{1}-1,n_{2}+1\}$ occur (Figure 3). The total number of distinct $\{n_{0},n_{1},n_{2}\}$ states can be enumerated via

[TABLE]

Due to the irreversibility of the binding process, at $t\to\infty$ we expect the system to consist only of $s_{\text{2}}\text{-network}$ strands: $P(n_{1},n_{2},t\to\infty)=0$ for all $\{n_{1},n_{2}\}\neq\{0,N_{\text{s}}\}$ and $P(0,N_{\text{s}},t\to\infty)=1$ as depicted in Figure 1g. We can obtain an alternate representation for Equation 11 by using the $n_{0}+n_{1}+n_{2}=N_{\rm s}$ constraint to represent $n_{2}$ so that $P(n_{0},n_{1},n_{2},t)\to P(n_{0},n_{1},t)$ and the master equation reads

[TABLE]

This representation is equivalent to Equation 11 and will be useful in deriving the distribution $P(n_{0},n_{1},t)$ from which $P(n_{1},n_{2},t)$ can be obtained. We now nondimensionalize our model by measuring time in units of the typical bond formation time, $\lambda^{-1}$. Henceforth, time $t$ will be dimensionless and $\lambda$ will no longer appear (equivalently, we set $\lambda=1$ in Equations 11 and 13). The mean number of strand types $\langle n_{\ell}(t)\rangle$ in a single micro-regions are defined by

[TABLE]

for $\ell=1,2$, under the $0\leq n_{1}+n_{2}\leq N_{\text{s}}$ constraint. The corresponding mass-action equations can be derived by multiplying Equation 11 by $n_{\ell}$ and by summing over $n_{1},n_{2}$ under the same constraint so that

[TABLE]

Equations 15a-c can be solved under the initial condition $n_{0}(0)=N_{\text{s}}$, representing all strands being unbound at $t=0$. We find

[TABLE]

so that $\langle n_{\ell}(t\to\infty)\rangle\to 0$ for $\ell=0,1$ and $\langle n_{2}(t\to\infty)\rangle\to N_{\text{s}}$. Equations 16a-c represent average values calculated across all micro-regions at time $t$ under quenched binding. We compare the $t\to\infty$ limit of Equations 16a-c to Equation 2 which estimates average strand numbers using combinatoric arguments. To do so, we evaluate $\langle m\rangle=\langle n_{1}\rangle+2\langle n_{2}\rangle$ to find

[TABLE]

from which we calculate the average extent of reaction $\langle p(t)\rangle={\langle m(t)\rangle}/{2N_{\text{s}}}$

[TABLE]

Inverting the transcendental Equations 17 and 18 for general $\alpha$ is not possible, however, under neutral cooperativity $\alpha=1$, we find

[TABLE]

A simple analysis of Equations 18 and 19 reveals that $\langle p(t)\rangle$ is a monotonically increasing function of $t$ for all $\alpha>0$, which is expected given that end-groups bind but do not unbind. For $\alpha=1$, Equations 16 can be recast as

[TABLE]

Equations 20 have the same form as Equation 2, obtained using mean-field arguments. This implies that the mean-field approach for a given extent of reaction $p$ and $\alpha=1$ corresponds to an irreversible (quenched) stochastic process halted at time $t^{*}$ such that $\langle p(t^{*})\rangle$ in Equation 19 satisfies $\langle p(t^{*})\rangle=1-e^{-t^{*}}=p$. We plot the normalized average strand numbers $\langle n_{\ell}(t)\rangle/N_{\text{s}}$ as a function of time and as derived from Equations 16a-c in Figure 4a-c, for $N_{\text{s}}=40$ and different values of $\alpha$. We find that $s_{\text{1}}$-strand formation is favored at smaller $\alpha$, and $s_{\text{2}}$-strand formation is favored at higher $\alpha$, as might be expected. Since $\langle p(t)\rangle$ is a monotonic function of time we can plot $\langle n_{\ell}(t)\rangle/N_{\text{s}}$ using Equations 16a-c as parametric equations against $\langle p(t)\rangle$ given in Equation 18. Results are shown in Figure 4d-f for various values of $\alpha$. These curves differ from those in Figure 2d obtained under equilibration and calculated via Equations 7a-c. Most noticeably, $\langle n_{1}(t)\rangle$ loses its symmetry about $\langle p(t)\rangle=0.5$ and becomes skewed.

The master Equation 11 also allows us to derive the time-dependent likelihood of each of the many possible configurations (enumerated in Equation 12), a much more powerful tool than average quantities. For example, Equation 11 can be solved to find $P(n_{1},n_{2},t)$ for all times $t$ once the initial condition is specified. We set this to be $P(n_{1}=0,n_{2}=0,t=0)=1$, and $P(n_{1},n_{2},t=0)=0$ for all other values of $n_{1},n_{2}\neq 0$, so that the micro-region is initially made only of free strands. If one chooses to solve Equation 13 to find $P(n_{0},n_{1},t)$ the equivalent initial conditions are $P(n_{0}=N_{\rm s},n_{1}=0,t=0)=1$ and $P(n_{0},n_{1},t=0)=0$ for $n_{0}\neq N_{\rm s}$ and $n_{1}\neq 0$. We solve Equation 13 for $P(n_{0},n_{1},t)$ rather than Equation 11 for $P(n_{1},n_{2},t)$ as the analytical computations are simpler. From $P(n_{0},n_{1},t)$ we can then derive $P(n_{1},n_{2},t)$ by changing variables through the $n_{0}+n_{1}+n_{2}=N_{\rm s}$ constraint. To proceed, we introduce the generating function $G(z_{0},z_{1},t)$ defined as

[TABLE]

under the constraint $0\leq n_{0}+n_{1}\leq N_{\text{s}}$. Upon multiplying Equation 13 by $z_{0}^{n_{0}}z_{1}^{n_{1}}$ and summing over $n_{0},n_{1}$, under the same constraint, we find the following differential equation for $G(z_{0},z_{1},t)$

[TABLE]

Equation 22 is coupled to the corresponding initial condition $G(z_{0},z_{1},t=0)=z_{0}^{N_{\text{s}}}$. Using the method of characteristics, we find

[TABLE]

After performing a Taylor series expansion in $z_{0},z_{1}$ and upon comparison with Equation 21 we find

[TABLE]

From the constraint $n_{0}=N_{\rm s}-n_{1}-n_{2}$ we can finally write

[TABLE]

Note that $P(n_{1},n_{2},t\to\infty)=0$ for $\{n_{1},n_{2}\}\neq\{0,N_{\text{s}}\}$ and that $P(0,N_{\text{s}},t\to\infty)=1$ as expected from a forward process. Figure 5 shows the probability of individual configurations $P(n_{1},n_{2},t)$ of micro-regions with $N_{\text{s}}=3$ plotted parametrically against $\langle p(t)\rangle$ for different values of $\alpha$. Two different configurations are highlighted: $\{n_{0},n_{1},n_{2}\}=$ $\{0,3,0\}$ and $\{1,0,2\}$. For $\alpha$ = 2, when end-group binding is cooperative, the probability of configurations with more $s_{\text{1}}\text{-strands}$ decreases and those with more $s_{\text{2}}\text{-strands}$ increases compared to the neutral ($\alpha=1$) or uncooperative ($\alpha=0.5$) cases shown here. In highly cooperative scenarios, once a network strand has bound the transition towards a fully-bound $s_{2}$-strand is fast. In Figure 6a-d we plot the micro-region configuration probabilities $P(n_{1},n_{2},t)$ as a function of the average extent of reaction $\langle p(t)\rangle$ with increasing $N_{\text{s}}$. The highest value of $N_{\text{s}}=40$ we used results in 861 distinct $\{n_{0},n_{1},n_{2}\}$ micro-region configurations, as per Equation 12, all with non-zero probability at finite time. Larger values of $N_{\text{s}}$ are possible, but graphically difficult to display.

By inserting the explicit expression for $P(n_{1},n_{2},t)$ from Equation 25 into Equation 14 we evaluate the average values $\langle n_{\ell}(t)\rangle$ for $\ell=1,2$ to reobtain the same expressions for $\langle n_{\ell}(t)\rangle$, for $\ell=0,1,2$ already displayed in Equations 16. From Equation 25 we can also calculate the second moments $\langle n^{2}_{\ell}(t)\rangle$ as

[TABLE]

for $\ell=1,2$, and the correlation function

[TABLE]

from which we can derive $\langle m^{2}(t)\rangle$

[TABLE]

Using Equations 16, 17, 26, 27 we find the variances

[TABLE]

for $\ell=1,2$. Finally, $\langle n^{2}_{0}\rangle$ is obtained as

[TABLE]

Equations 30 and 16 yield \mathrm{Var}\big{(}n_{0}(t)\big{)}=\langle n^{2}_{0}(t)\rangle-\langle n_{0}(t)\rangle^{2}. Explicit expressions for $\langle n^{2}_{\ell}(t)\rangle$, $\langle m^{2}(t)\rangle$, \mathrm{Var}\big{(}n_{\ell}(t)\big{)} and \mathrm{Var}\big{(}m(t)\big{)} are presented in Section A of the Appendix. In Figure 6e-h we show the parametric plots of the average strand fractions $\langle n_{\ell}(t)\rangle/N_{\rm s}$ for $\ell=0,1,2$ against the average extent of reaction $\langle p(t)\rangle$ for several values of $N_{\rm s}$. The associated standard deviations calculated as the square root of the variance in Equation 29 are also displayed. As can be expected, fluctuations decrease as $N_{\rm s}$ increases.

Figure 7a shows the parametric plot of \mathrm{Var}\big{(}n_{2}(t)\big{)} against $\langle p(t)\rangle$ for different values of $\alpha$ and $N_{\text{s}}=40$. For strong uncooperative binding ($\alpha\to 0$) and small extents of reaction $\langle p(t)\rangle$, only free strands bind to the network and the variance is very small. However, once all strands have bound at least at one end, and $\langle p(t)\rangle\sim 0.5$, the dangling strands transition to the fully bound state and the variance starts increasing. In Equation 41 of Section A of the Appendix we give an exact analytical expression for \mathrm{Var}\big{(}n_{2}(t)\big{)}; a simple calculation shows that the maximum variance is $N_{\rm s}/4$ for all values of $\alpha$, and is attained at smaller average extents of reaction $\langle p(t)\rangle$ as $\alpha$ increases. In Figure 7b we plot the parametric dependence of \mathrm{Var}\big{(}p(t)\big{)} against the average extent of reaction $\langle p(t)\rangle$ with variable $\alpha$ and $N_{\text{s}}=40.$ For very small values of $\alpha\to 0$, \mathrm{Var}\big{(}p(t)\big{)} is bimodal and approximately zero at $\langle p(t)\rangle\to 0.5$. This is because, as discussed above, for $\alpha\to 0$ end-group binding occurs only on free strands for $\langle p(t)\rangle<0.5$ and the most likely configurations are those with $s_{0}$ and $s_{1}$ strands. As $\langle p(t)\rangle\to 0.5$, only $s_{1}$ dangling strands remain so that $\langle n_{1}(t)\rangle\to N_{\rm s}$, $\langle n^{2}_{1}(t)\rangle\to N^{2}_{\rm s}$, \mathrm{Var}\big{(}n_{2}(t)\big{)}\to 0 and as a result, \mathrm{Var}\big{(}p(t)\big{)}\to 0. Fully bound strands start emerging for $\langle p(t)\rangle>0.5$, increasing \mathrm{Var}\big{(}p(t)\big{)}. As $\alpha$ increases the variance increases for all $\langle p(t)\rangle$ and the minimum at $\langle p(t)\rangle\sim 0.5$ turns into a maximum. A more detailed discussion is presented in Section B of the Appendix.

In Figure 7c we plot \mathrm{Var}\big{(}n_{2}(t)\big{)}/N^{2}_{\rm s} against $\langle p(t)\rangle$ for different values of $N_{\text{s}}$; the curves decrease in magnitude as $N_{\rm s}$ increases. Finally, in Figure 7d we plot $\langle n_{1}(t)\rangle/N_{\rm s}$ parametrically against $\langle p(t)\rangle$ for several values of $\alpha$. The curves decrease with $\alpha$ once $\langle p(t)\rangle$ is fixed. This also follows from Equation 16b which implies that $\langle n_{1}(t)\rangle$ is a decreasing function of $\alpha$ for any time $t$. Since $\langle p(t)\rangle$ is univocally associated to $t$ via Equation 18 it also follows that $\langle n_{1}(t)\rangle$ is a decreasing function of $\alpha$ for any $\langle p(t)\rangle$. Equations 16b and 18 reveal that the maximum $\langle n_{1}(t_{\rm max})\rangle/N_{\rm s}=(\alpha/2)^{(1-2/\alpha)}$ is attained at $\langle p(t_{\rm max})\rangle=1-((\alpha+1)/\alpha)e^{-2t_{\rm max}}$ where $t_{\rm max}=\ln(2/\alpha)/(2-\alpha)$. One can easily verify that $\langle p(t_{\rm max})\rangle$ is a decreasing function of $\alpha$ as well. These results can be expected as larger $\alpha$ favors the formation of fully bound strands. Thus, for a given average extent of reaction $\langle p(t)\rangle$ the fraction of dangling ends decreases with $\alpha$, and the maximum is found at an average extent of reaction $\langle p(t)\rangle$ that also decreases with $\alpha$.

II.2.2 Dynamic end-group rearrangement/redistribution

We now consider an equilibration process that allows the bound end-groups in a micro-region to dynamically rearrange, attaching and detaching until thermodynamic equilibrium is reachedBowman and Kloxin (2012) while maintaining a fixed total number of $m$ bound end-groups. We assume that $m<2N_{\text{s}},(p<1)$ so that the reaction is not complete and multiple $\{n_{0},n_{1},n_{2}\}$ configurations are possible. Experimental realizations include the formation reversible hydrazone bonds Roberts et al. (2007), imine bonds Yang et al. (2012), or guest-host interactions Rodell et al. (2015). Here, an intact, $s_{\text{2}}\text{-}$strand may detach at one of its ends to form a dangling $s_{\text{1}}\text{-}$strand, while a free $s_{\text{0}}\text{-}$strand binds to the network to form another $s_{\text{1}}\text{-}$strand. The reverse process where two $s_{\text{1}}\text{-}$strands become an $s_{\text{2}}\text{-}$ and $s_{\text{0}}\text{-}$strand is also possible. In all scenarios $m=n_{1}+2n_{2}$ is fixed, but there are many distinct ways for the bound end-groups to distribute in $s_{\text{1}}\text{-}$ or $s_{\text{2}}\text{-}$ strands. The final equilibrium configuration is independent of initial conditions so our results will depend only on the selected value of $m$. We write the reversible master equation for $P(n_{1},n_{2},t)$ as

[TABLE]

Here, $\kappa$ is dimensionless and represents the rearrangement rate measured in terms of the binding rate. The first term on the right hand side of Equation 31 accounts for the formation of a fully bound $s_{\text{2}}\text{-strand}$ and a free $s_{\text{0}}\text{-strand}$ from two dangling $s_{\text{1}}\text{-strands}$ ($2s_{\text{1}}\rightarrow s_{\text{0}}+s_{\text{2}}$). Here, the bound end-group of one of the two $s_{\text{1}}\text{-strands}$ exchanges with the unbound end-group of the other leading to the $\{n_{0}-1,n_{1}+2,n_{2}-1\}\to\{n_{0},n_{1},n_{2}\}$ transition. The combinatorial factor enumerates the number of $s_{\text{1}}\text{-strands}$ present in the micro-region and the 2 prefactor represents both $s_{\text{1}}\text{-strands}$ being able to exchange with the other. Finally, the reactivity parameter $\alpha$ is squared, since two dangling ends must bind to form a fully bound strand. The second term represents a fully bound $s_{\text{2}}\text{-strand}$ detaching on one end while promoting the binding of a free $s_{\text{0}}\text{-strand}$, giving rise to two dangling $s_{\text{1}}\text{-strands}$. This process is represented by the $\{n_{0}+1,n_{1}-2,n_{2}+1\}\to\{n_{0},n_{1},n_{2}\}$ transition. The factors $(N_{\text{s}}-n_{1}-n_{2}+1)$ and $(n_{2}+1)$ represent the number of $s_{\text{0}}\text{-}$ and $s_{\text{2}}\text{-}$network strands available, respectively. The prefactor 4 accounts for the number of possible bond movements: either of the two bound end-groups on the $s_{\text{2}}$-strand can relocate to either of the two unbound end-groups of the $s_{\text{0}}$-strand, yielding a total of four combinations. The last two terms represent the same two processes described above, but driving the system away from $\{n_{0},n_{1},n_{2}\}$. Note that there are no terms in Equation 31 representing bonds leaving an $s_{\text{2}}$-strand to populate an $s_{\text{1}}$-strand; this transition would not change the overall the micro-region configuration $\{n_{0},n_{1},n_{2}\}$. The probability $P_{\text{b}}(m,t)$ of having $m$ bound-ends at time $t$ can be written as

[TABLE]

where all possible $n_{1},n_{2}$ combinations that yield $m=n_{1}+2n_{2}$ are included. Using Equation 31 it can be easily verified that $P_{\text{b}}(m,t)=P_{\rm b}(m,t=0)$. As expected, the master Equation 31 only rearranges the distribution of $s_{\text{1}}$ and $s_{\text{2}}$ strands but $m$ remains unchanged. We thus assume the system is initiated with a given $m$ so that $n_{1}+2n_{2}=m$ at all times. In addition to this constraint, the number of strands is fixed so that $n_{0}+n_{1}+n_{2}=N_{\text{s}}$. We can thus cast Equation 31 in terms of only one of the $n_{0},n_{1}$ or $n_{2}$ populations. We choose $n_{2}$ and determine the steady state $P(n_{1},n_{2},t\to\infty)\equiv P^{*}(n_{2})$ by imposing detailed balance between the first and the last term on the right hand side of Equation 31, or equivalently, the second and the third, since it can be easily verified that the conditions are the same. We find

[TABLE]

which can be solved to yield

[TABLE]

where $Z_{m,N_{\text{s}}}$ is the normalization constant

[TABLE]

This result is the same as Equation 8: the combinatoric approach for a fixed $p$ is equivalent to allowing for relaxation on the network with a fixed number of bound ends and $m=2pN_{\rm s}$.

II.2.3 End-group rearrangement/redistribution and bond formation

We now consider the two processes of bond formation and redistribution occurring simultaneously and combine the two master Equations 11 and 34 so that

[TABLE]

Fast annealing, fast binding, and quenched/irreversible binding are modeled by setting $\kappa\gg 1$, $\kappa\ll 1$, and $\kappa=0$, respectively.

Although a full analytical time-dependent solution can not be found, the effects of annealing can be observed in Figure 8. Here we parametrically plot $P(n_{1},n_{2},t)$ against $\langle p(t)\rangle$ using Equation 36 for $N_{\rm s}=m=5$, $\alpha=0.5$ and $\alpha=2$, under fast annealing ($\kappa=1000$) or quenched binding ($\kappa=0$). Since the rearrangement process allows for more configurations to be explored we expect cooperative effects to be more pronounced under fast annealing, than under quenched binding. In Figure 8a we set $\alpha=0.5$; since binding is uncooperative, annealing favors configurations with lower values of $n_{2}$. Indeed, the $k=1000$ curves show an increase in $P(5,0,t)$ compared to the corresponding $\kappa=0$ curves, whereas $P(3,1,t),P(1,2,t)$ decrease. Similar trends are observed in Figure 8b, where we set $\alpha=2$. Cooperative binding increases the likelihood of configurations with higher $n_{2}$, so in this case $P(1,2,t)$ increases while $P(3,1,t),P(5,0,t)$ decrease. Note that $P(5,0,t)$, $P(3,1,t)$, and $P(1,2,t)$ all obey the constraint $n_{1}+2n_{2}=N_{\rm s}=5$. For $\alpha=2$, and under quenched binding at $\kappa=0$, Equation 25 yields $P(5,0,t)=32e^{-10t}t^{5}$ which is maximized at $t=1/2$ corresponding to $\langle p(t=1/2)\rangle=1-3/2e\neq 1/2$, as per Equation 18. Similarly, $P(3,1,t)=8e^{-8t}t^{3}(1-(2t+1)e^{-2t})$ and $P(1,2,t)=2e^{-6t}t(1-(2t+1)e^{-2t})$ are also maximized at times that correspond to $\langle p(t)\rangle\neq 1/2$. None of the three distribution curves are thus symmetric about $\langle p(t)\rangle=1/2$. When $\kappa=1000$ however the master Equation 36 yields numerical results that are closely aligned with those derived from Equation 31 upon setting $n_{1}+2n_{2}=N_{\rm s}=5$. This is because annealing is much faster than binding and the time between binding events is much longer than the time for equilibration of a fixed number of bound strands. As a result, once a strand binds, the network can almost fully equilibrate before the next binding event occurs. The curves in Figure 8b for $\kappa=1000$ thus mirror Equation 34, with the proportions $P(5,0)$:$P(3,1)$:$P(1,2)$ following Equation 33 and become symmetric about $\langle p(t)\rangle=0.5$ as predicted by Equation 34 when $m=N_{\rm s}$. The same trends arise when comparing the quenched binding and the fast annealing curves for the uncooperative ($\alpha=0.5$) case.

In Figure 9 we plot $\langle n_{\ell}(t)\rangle/N_{\text{s}}$ parametrically against $\langle p(t)\rangle$ for $\alpha=$ 0.5, 1, 2 using the probability distribution in Equation 36 for $\ell=0,1,2$ and $\kappa=1000$. In all cases, the solutions closely match those obtained from the equilibrated distribution in Equation 34 for all values of $\alpha$ as can be seen upon comparison with Figure 2d. The most notable feature is the symmetry of $\langle n_{1}(t)\rangle$ about $\langle p(t)\rangle=0.5$ for all $\alpha$, a feature of the combinatoric approach, as discussed in Section II.1, Intermediate values of $\kappa\approx 1$ yield curves that interpolate between the two extremes $\kappa\gg 1$ and $\kappa=0$ shown here. Our results imply that networks formed via quenched end-group binding, as per Equation 25, should not be described by models that assume network strands equilibration via redistribution, as per Equations 1 and 9.

III Applications and Conclusions

In this work we studied the stochastic properties of bifunctional network strands that undergo an end-linking gelation process. We developed and analyzed a master equation to describe quenched and annealed binding events in micro-regions within a larger polymer network, and include a reactivity parameter to model cooperative effects. While typical models quantify “average” mean-field properties, we are able to evaluate the full probability distribution for any given configuration as a function of time and extent of reaction. By modeling the probability of a configuration within a micro-region, we can propose a crude framework to describe the effects of heterogeneity across the entire sample.

For example, our approach can be used in several polymer network applications under the assumption that a macroscopic region is comprised of a collection of statistically identical, independent, smaller micro-regions. For example, nano/micrometer scale differences in the polymer network properties can affect the fate of cells that are cultured on them Yang et al. (2016) as well as the mechanical properties of high-performance materials Li and Strachan (2016). If these properties depend on the local number of free, dangling, and intact strands, we can use the relevant probability distribution to evaluate the likelihood of a given configuration $\{n_{0},n_{1},n_{2}\}$ in any number of micro-regions sampled by e.g., a cell. The statistical distribution for each micro-region can then be used to construct the probability distribution of the entire macro-system and thus to estimate the chemical or mechanical properties of the polymer network, including their local variability.

Similar considerations can be applied to the study of elastic properties, in particular within phantom network theory which posits that the shear modulus of an ideal network depends on the number density of elastically effective network strands. Our results are readily applicable if we assume that all $s_{\text{2}}\text{-strands}$ in our models are elastically effective and the number of branchpoints is fixed. Starting from the probability $P(n_{1},n_{2},t)$ that a micro-region is in the $\{n_{0},n_{1},n_{2}\}$ configuration, we can also compute the likelihood that a given threshold is met, say, $n_{2}\geq n_{2}^{*}$. This quantity can then be interpreted as the probability for a “bond” to stretch across a micro-region. One can then calculate the likelihood that a given number of contiguous micro-regions with $n_{2}\geq n_{2}^{*}$ span the sample through percolation, leading to a dramatic stiffening of the network.

Finally, our work can also be applied to the study of network degradation, which has attracted much attention over the past two decades as degradable sites have been increasingly incorporated into experimental realizations of end-linking polymerization. These degradable strands typically cleave by enzymatic, hydrolytic, photolytic, or other chemical mechanisms and allow for a reverse gelation process. Here, strands initially exist in the fully bound, intact state where both ends are unreacted. Reverse gelation occurs via reaction or degradation of either end, so that intact strands first become dangling strands, and dangling strands then become free strands. Halting the extent of reaction is common in degradable networks as a way to tune the gel mechanics and this results in a large variability of the micro-region composition. Photodegradable networks Griffin and Kasko (2012); Norris et al. (2017); Xue et al. (2014); Käpylä et al. (2016); Norris et al. (2019); Wong et al. (2010), where end-groups are degraded by exposure to light, are of particular interest as they are uniquely suited to spatially pattern network stiffness, with a high degree of control Norris et al. (2016). Some mathematical models of reverse gelation have been formulated by adapting models of gelation Metters et al. (2000); more specific mean-field photodegradable network models have also been proposed Tibbitt et al. (2013); Norris et al. (2017). The present work can be adapted to model degradable networks by associating intact network strands to $s_{\text{0}}\text{-strands}$ (0 degraded end-groups), dangling strands to $s_{\text{1}}\text{-strands}$ (1 degraded end-group), and free strands to $s_{\text{2}}\text{-strands}$ (2 degraded end-groups). Cooperative effects arise in this context as the un-degraded end-groups of an intact strand might more readily react due to tension across the strand induced by network swelling. Once one of the end-groups has cleaved, and the strand dangles, the stress is removed so that the remaining un-degraded end-group is less susceptible to further degradation. Using our stochastic framework, one can calculate the probability of any given micro-region configuration, distinguishing between quenched and annealed network de-gelation reactions. Melting and collapse of rigidity can be then be described using percolation concepts.

Acknowledgements.

Funding for this work was provided by the National Institutes of Health through the NIH Director’s New Innovator Award Program, 1-DP2-OD008533. S.C.P.N gratefully acknowledges support from a Ruth L. Kirschstein Predoctoral Fellowship (NIH-F31DE026356). T. C. acknowledges support from the NSF through grant DMS-1814364. M.R.D. acknowledges support from the NSF through grant DMS-1814090 and Army Research Office W911NF-16-1-0165 . Both T.C. and M.R.D. also acknowledge support from the Army Research Office (W911NF-18-1-0345).

Appendix A Second moments

We here derive $\langle n^{2}_{\ell}(t)\rangle$, $\langle m^{2}(t)\rangle$, \mathrm{Var}\big{(}n_{\ell}(t)\big{)} and \mathrm{Var}\big{(}m(t)\big{)} for $\ell=0,1,2$ using the explicit form for $P(n_{1},n_{2},t)$ as given in Equation 25. We begin with

[TABLE]

and the associated variances for $\ell=1,2$. Using the binomial theorem we find

[TABLE]

which, coupled with Equation 16b for $\langle n_{1}(t)\rangle$ leads to

[TABLE]

Similarly, Equation 37 for $\ell=2$ yields

[TABLE]

which coupled with Equation 16c for $\langle n_{2}(t)\rangle$ leads to

[TABLE]

Equation 41 is maximized for time $t_{\rm M}$ implicitly given by

[TABLE]

which corresponds to \mathrm{Var}\big{(}n_{2}(t_{\rm M})\big{)}=N_{\rm s}/4, independent of the value of $\alpha$. To evaluate $\langle m^{2}(t)\rangle$ we must first calculate the correlation function

[TABLE]

which yields

[TABLE]

We can now evaluate $\langle m^{2}(t)\rangle$ using Equation 28, 38 and A

[TABLE]

from which the variance is obtained as

[TABLE]

where we used Equations 16b-c to evaluate $\langle m(t)\rangle$. Finally, using the constraint $n_{0}=N_{\rm s}-n_{1}-n_{2}$ we find

[TABLE]

which together with Equation 16a finally yields

[TABLE]

Appendix B Strong uncooperative binding

All evaluations in the main text assume $\alpha\neq 0$, since the completely uncooperative case ($\alpha=0$) would not allow for the formation of $s_{2}$ strands. Setting $\alpha=0$ however can be used to explore the short time behavior when $\alpha\to 0$. This is the case of highly uncooperative binding where although rare, the formation of a fully bound $s_{2}$ strand is still possible. Setting $\alpha=0$ in Equations 16b-c and 18, so that $e^{-\alpha t}\to 1$ for all times, we find

[TABLE]

Upon setting $\alpha=0$ in Equation 46 we also find

[TABLE]

Note that Equation 51 yields $\langle p(t)\rangle<0.5$ for all times, implying that for $\alpha=0$ the reaction cannot be completed, as expected since fully bound strands cannot emerge. Equation 52 also reveals that \mathrm{Var}\big{(}p(t)\big{)} is symmetric about $\langle p\rangle=1/4$ and its maximum is attained at \mathrm{Var}\big{(}p(t)\big{)}=(16N_{\rm s})^{-1}. The above results still apply in the $\alpha\to 0$ limit, albeit for $\alpha t\ll 1$ where $e^{-\alpha t}\to 1$. For example \mathrm{Var}\big{(}p(t)\big{)} follows Equation 52 in Figure 7b up to $\langle p(t)\rangle\sim 0.5$. At longer times, since $\alpha$ is small but not zero, the binding will proceed, and $s_{2}$ strands will emerge. We can thus reevaluate Equations 16b-c and 18, for $\alpha\to 0$ but at long times where $e^{-2t}\to 0$ and $e^{-\alpha t}\neq 0$ so that

[TABLE]

Finally, in the $e^{-2t}\to 0$ limit, Equation 46 becomes

[TABLE]

The results in Equations 55 and 56 indicate that as $t\to\infty$, $\langle p(t)\rangle>0.5$ and $\langle p(t)\rangle\to 1$. Furthermore, we observe that the shape of \mathrm{Var}\big{(}p(t)\big{)} in Equation 56 is the same as in Equation 52 as also emerges from the bimodal form in Figure 7b. Finally, we note that for $t\to\infty$, $\langle n_{1}(t)\rangle\to 0$, even as $\alpha\to 0$ since eventually all strands will be fully bound and $\langle n_{2}(t)\rangle\to N_{\rm s}$.