Multistate reversible copolymerization of nonMarkovian chains under low conversion conditions

TL;DR

This paper analytically solves the reversible multistate copolymerization kinetics under low conversion conditions, revealing how non-Markovian chain sequences are characterized by matrices linked to monomer units, affecting growth and thermodynamics.

Contribution

It provides an analytical solution for the kinetics of multistate reversible copolymerization, introducing a matrix-based characterization of non-Markovian sequence properties.

Findings

Matrices determine growth velocity and sequence statistics.

Non-Markovian sequences are characterized by state-dependent matrices.

Thermodynamics of copolymerization are derived from the matrices.

Abstract

The reversible kinetics of copolymerization is solved analytically for the multistate mechanism proposed by B. D. Coleman and T. G. Fox [J. Chem. Phys. 38, 1065 (1963)] under low conversion conditions where the concentrations of monomeric species are chemostatted and stay constant in time. Although the rates of this mechanism only depends on the currently attached or detached monomer, the growing macromolecular chain forms a nonMarkovian sequence that is characterized by matrices associated with every monomeric unit composing the sequence. These matrices are obtained by solving the kinetic equations and they determine the growth velocity of the copolymers, the statistical properties of its possible sequences, as well as the thermodynamics of the copolymerization process.

Tables

Table 1. Table 1: Growth of a copolymer in the conditions ( 108 ) with c A = 5 subscript 𝑐 A 5 c_{\rm A}=5 : Comparison between the values of the tip probabilities ( 46 ) obtained by theory with 10 4 superscript 10 4 10^{4} iterations of Eq. ( 42 ) and numerically by Gillespie’s algorithm generating 10 7 superscript 10 7 10^{7} copolymer sequences, each of total time t = 100 𝑡 100 t=100 . Under the chosen conditions, the mean growth velocity is v = 14.580 𝑣 14.580 v=14.580 , the free-energy driving force ε = 1.8623 𝜀 1.8623 \varepsilon=1.8623 , the sequence disorder D = 0.4534 𝐷 0.4534 D=0.4534 , and the affinity A = ε + D = 2.3157 𝐴 𝜀 𝐷 2.3157 A=\varepsilon+D=2.3157 . The mean length of the sequences is thus equal to ⟨ l ⟩ t = v ​ t = 1458 subscript delimited-⟨⟩ 𝑙 𝑡 𝑣 𝑡 1458 \langle l\rangle_{t}=vt=1458 . The same simulation data are used for Fig.  1 . The numbers are rounded off to five decimal digits.

$m_{l-2}$	$m_{l-1}$	$m_l$	$i$	theory	simulation
			$1$	$0.66667$	$0.66671$
			$2$	$0.33333$	$0.33329$
		A	$1$	$0.52263$	$0.52262$
		B	$1$	$0.14403$	$0.14409$
		A	$2$	$0.30048$	$0.30048$
		B	$2$	$0.03285$	$0.03281$
	A	A	$1$	$0.41483$	$0.41494$
	B	A	$1$	$0.10781$	$0.10768$
	A	B	$1$	$0.11328$	$0.11331$
	B	B	$1$	$0.03076$	$0.03078$
	A	A	$2$	$0.26832$	$0.26832$
	B	A	$2$	$0.03217$	$0.03216$
	A	B	$2$	$0.02889$	$0.02888$
	B	B	$2$	$0.00396$	$0.00393$
A	A	A	$1$	$0.33315$	$0.33337$
B	A	A	$1$	$0.08168$	$0.08158$
A	B	A	$1$	$0.08518$	$0.08512$
B	B	A	$1$	$0.02262$	$0.02257$
A	A	B	$1$	$0.09018$	$0.09020$
B	A	B	$1$	$0.02309$	$0.02311$
A	B	B	$1$	$0.02422$	$0.02422$
B	B	B	$1$	$0.00654$	$0.00655$
A	A	A	$2$	$0.23791$	$0.23788$
B	A	A	$2$	$0.03041$	$0.03044$
A	B	A	$2$	$0.02785$	$0.02786$
B	B	A	$2$	$0.00431$	$0.00430$
A	A	B	$2$	$0.02532$	$0.02534$
B	A	B	$2$	$0.00357$	$0.00355$
A	B	B	$2$	$0.00336$	$0.00333$
B	B	B	$2$	$0.00060$	$0.00059$

Equations

$$\displaystyle m_{1}m_{2}\cdots m_{l-1}^{*}\ +\ m_{l}\ \$$

$$\displaystyle m_{1}m_{2}\cdots m_{l-1}^{**}\ +\ m_{l}\ \$$

$$\displaystyle m_{1}m_{2}\cdots m_{l-1}m_{l}^{*}\ \$$

$$\frac{d}{dt}{P}_t(m_1\cdots m_l,i)=w_{i,+m_l}{P}_t(m_1\cdots m_{l-1},i)$$

$$+\sum _{m_{l+1}=1}^M{w}_{i,-m_{l+1}}{P}_t(m_1\cdots m_l{m}_{l+1},i)$$

$$+w_{j\to i}{P}_t(m_1\cdots m_l,j)$$

$$\displaystyle-\bigg{(}w_{i,-m_{l}}+\sum_{m_{l+1}=1}^{M}w_{i,+m_{l+1}}+w_{i\to j}\bigg{)}P_{t}(m_{1}\cdots m_{l},i),$$

$$p_t(i,l)$$

$$p_t(m_l,i,l)$$

$$p_t(m_{l-1}{m}_l,i,l)$$

$$a_i\equiv \sum _m{w}_{i,+m},$$

$$\frac{d}{dt}{p}_t(i,l)=a_i{p}_t(i,l-1)$$

$$+\sum _{m_{l+1}}{w}_{i,-m_{l+1}}{p}_t(m_{l+1},i,l+1)+w_{j\to i}{p}_t(j,l)$$

$$-\sum _{m_l}{w}_{i,-m_l}{p}_t(m_l,i,l)-(a_i+w_{i\to j})p_t(i,l),$$

$$\frac{d}{dt}{p}_t(m_l,i,l)=w_{i,+m_l}{p}_t(i,l-1)$$

$$+\sum _{m_{l+1}}{w}_{i,-m_{l+1}}{p}_t(m_l{m}_{l+1},i,l+1)+w_{j\to i}{p}_t(m_l,j,l)$$

$$-(w_{i,-m_l}+a_i+w_{i\to j})p_t(m_l,i,l),$$

$$⋮$$

$$p_t(m_{l-r+1}\cdots m_{l-1}{m}_l,i,l)$$

$$=\exp (s_{q}t+\text{i}ql){\chi }_q(m_{l-r+1}\cdots m_{l-1}{m}_l,i)$$

$$s_q=-\text{i}qv-\mathcal{D}{q}^2+O(q^3),$$

$$\psi (m_{l-r+1}\cdots m_{l-1}{m}_l,i)={\chi }_0(m_{l-r+1}\cdots m_{l-1}{m}_l,i).$$

$$0=-w_{i\to j}\psi (i)+w_{j\to i}\psi (j),$$

$$0=w_{i,+m_l}\psi (i)+\sum _{m_{l+1}}{w}_{i,-m_{l+1}}\psi (m_l{m}_{l+1},i)$$

$$-(w_{i,-m_l}+a_i+w_{i\to j})\psi (m_l,i)+w_{j\to i}\psi (m_l,j),$$

$$⋮$$

$$-\text{i}v\psi (i)$$

$$\psi (1)+\psi (2)=1.$$

$$v$$

$${\bm{\mathsf{W}}}_{0}\equiv\left(\begin{array}[]{cc}-w_{1\to 2}&w_{2\to 1}\\ w_{1\to 2}&-w_{2\to 1}\end{array}\right),$$

$${\bm{\mathsf{W}}}_{\pm m}\equiv\left(\begin{array}[]{cc}w_{1,\pm m}&0\\ 0&w_{2,\pm m}\end{array}\right),$$

$$\mathsf{A}\equiv \sum _m{\mathsf{W}}_{+m}.$$

$$\mathbf{\Psi }$$

$$\mathbf{\Psi }(m_l)$$

$$\mathbf{\Psi }(m_{l-1}{m}_l)$$

$$\sum _{m_l}\mathbf{\Psi }(m_l)$$

$$\sum _{m_{l-1}}\mathbf{\Psi }(m_{l-1}{m}_l)$$

$$0$$

$$\psi (1)=\frac{w_{2\to 1}}{w_{1\to 2}+w_{2\to 1}},$$

Full text

I Introduction

In 1963, Coleman and Fox proposed a kinetic mechanism for the growth of polymers with nonMarkovian diastereosequence distributions.CF63JCP ; CF63JACS Their work concerns the tacticity of polymers, i.e., the property that consecutive identical monomeric units composing a polymer may have isotactic or syndiotactic placements, possibly forming atactic sequences.CF63JPS In the proposed mechanism, the growing polymer has two possible reactive states, each with its own stereospecificity for the attachment of monomers. Moreover, the attachment rates are assumed to be independent of previously incorporated monomers.CF63JCP If the polymer always remained in a single reactive state during the growth, the sequence would form a Bernoulli chain. Therefore, the possibility for the reactive state of the growing polymer to randomly change between two or more reactive states modifies the statistical properties of the chain and allows correlations to appear between consecutive monomeric units, generating instead a nonMarkovian chain. In the fully irreversible growth regime where the detachment rates are negligible in front of the attachment rates, Coleman and Fox showed that the nonMarkovian chains generated by this multistate mechanism have sequence probability distributions given by products of matrices.CF63JCP

Now, the multistate mechanism may also be considered for coplymerization processes, in which different species of monomers are incorporated in the growing chain. Moreover, in order to investigate the thermodynamic properties of such processes, the reversed reactions of monomeric detachments should no longer be neglected, as previously studied by numerical simulations in Ref. AG09, . In this context, the question arises whether the analytical methods developed by Coleman and FoxCF63JCP can be extended to reversible multistate processes.

The purpose of the present paper is to show that such an extension can indeed be achieved and that the sequence probability distributions are also given by matrix products in the general case of reversible multistate copolymerization processes under low conversion conditions.CF63JCP These latter are required in order to achieve the stationarity of monomer concentrations during the copolymer growth. If stationarity did not hold, the copolymer sequence would depend on the peculiar time variation undergone by the monomer concentrations in the solution,KZMI18 a situation encountered for instance in closed reactors where the population of monomers is depleted as the copolymers are growing until full equilibrium is reached.SP85 ; S87 ; S92 ; BL17 Instead, it is here supposed with Coleman and FoxCF63JCP that the monomer concentrations stay constant in time, which is the case under low conversion conditions or in an open reactor where the solution is continuously fed with monomers. The process is described using a stochastic approach at the level of a single copolymer. In this approach, the time evolution is considered for the probability to find the copolymer with given sequence, length, and reactive state. The time dependence of this probability is ruled by kinetic equations in terms of the transition rates associated with the different reactions in the process.

The content of the paper is the following.

In Section II, the kinetic equations of multistate reversible copolymerization are introduced and solved in the long-time limit in terms of matrices of size equal to the number of reactive states in the mechanism. A matrix is associated with every species of monomeric units in the growing copolymer sequence, as well as with the stationary probability distribution of the different reactive states of the copolymer. The probability of any given sequence is thus obtained as a corresponding matrix product. The mean growth velocity is shown to be given in terms of a velocity matrix, satisfying a self-consistent equation, which plays a central role in the theory. Furthermore, these different matrices determine the probability of any given monomeric subsequence in the bulk of the whole sequence, as well as the mean value, the variance, and the correlation function of any observable quantity defined along the sequence. The results previously obtained in the fully irreversible regime by Coleman and FoxCF63JCP are recovered.

In Section III, thermodynamics is developed for multistate reversible copolymerization in steady growth regimes.AG08 In particular, the expression of the entropy production rate is deduced and the consequences of the equilibrium detailed balance conditions are studied.

In Section IV, theory is compared with simulations using Gillespie’s algorithmG76 ; G77 for two illustrative examples: on the one hand, an example where the thermodynamic equilibrium limit exists and, on the other hand, one without equilibrium.

Conclusion and perspectives are given in Section V.

II Kinetics

II.1 The kinetic equations

We consider the Coleman-Fox multistate mechanismCF63JCP including the reversed reactions for the copolymerization of chains composed of different species of monomeric units $m=1,2,...,M$:

[TABLE]

in the case of two reactive states, which are here below denoted by $i=1,2$. The transition rates of the reactions (1) and (2) are $w_{i,\pm m_{l}}$ and those of the transitions between the two states are $w_{1\to 2}$ and $w_{2\to 1}$. The attachment rates are proportional to the concentration $c_{m}$ of the corresponding monomer in the solution, $w_{i,+m}=k_{i,+m}c_{m}$, with some rate constant $k_{i,+m}$. Instead, the detachment rates, as well as the rates of transitions between the two states, do not depend on the monomer concentrations, $w_{i,-m}=k_{i,-m}$ and $w_{i\to j}=k_{ij}$. Under low conversion conditions, we may assume that the monomer concentrations $c_{m}$ and thus the rates stay constant in time.CF63JCP

The copolymerization of a single chain is described as a stochastic process in terms of the probabilities $P_{t}(m_{1}m_{2}\cdots m_{l},i)$ to find the polymer with the different possible sequences $m_{1}m_{2}\cdots m_{l}$ of length $l$ and states $i$ at the time $t$. These probabilities are ruled by the following infinite hierarchy of coupled kinetic equations,

[TABLE]

where $j=2$ if $i=1$ or $j=1$ if $i=2$. The total probability is conserved in time, $\sum_{\omega}P_{t}(\omega)=1$, where the sum extends over the sequences $\omega=m_{1}m_{2}\cdots m_{l}$ and the states $i$.

II.2 Solving the kinetic equations

The kinetic equations are solved by introducing the following set of probabilities:

[TABLE]

Introducing the notation

[TABLE]

these probabilities obey the following equations as a consequence of the kinetic equations (4):

[TABLE]

and equations similar to Eq. (10) for the probabilities (II.2) and the next ones. We note that we recover Eq. (9) by summing Eq. (10) over $m_{l}$, and similarly for the further equations in the hierarchy.

Since these equations are linear, their general solution can be written as a linear superposition of solutions of the form

[TABLE]

with an arbitrary parameter $-\pi<q\leq+\pi$. The exponential rate $s_{q}$ is here expected to have the following form,

[TABLE]

where $\i=\sqrt{-1}$, $v$ is the mean growth velocity of the copolymer chain counted in monomers per second, and $\cal D$ is the diffusivity of the random drift of the length $l$. The particular solutions (11) are substituted in Eqs. (9), (10),… of the hierarchy. We carry out an expansion in powers of $q$ around $q=0$ and use the fact that $s_{0}=0$ together with the notation

[TABLE]

At order $q^{0}$, the following equations are obtained

[TABLE]

and equations similar to Eq. (15) for the quantities (13) with the next values of $r$.

At order $q^{1}$, Eq. (9) give

[TABLE]

where $j=2$ if $i=1$ or $j=1$ if $i=2$, $m=m_{l+1}$, and $\psi^{\prime}(i)=d\chi_{q}(i)/dq|_{q=0}$. Moreover, we suppose that the following normalization condition is satisfied:

[TABLE]

Summing Eq. (16) over $i=1,2$, we obtain the relation giving the mean growth velocity

[TABLE]

In order to solve the previous equations, we introduce a $2\times 2$ matrix describing the coupling between the pairs of equations due to the rates $w_{i\to j}$ of the transitions between the reactive states

[TABLE]

as well as other $2\times 2$ matrices with the rates of the attachment and detachment reactions

[TABLE]

together with

[TABLE]

Furthermore, we also define the $2\times 2$ matrices

[TABLE]

As before, we have that

[TABLE]

With these definitions, Eqs. (14), (15),… of the hierarchy can be rewritten in the following matricial form:

[TABLE]

First, we note that Eqs. (14) and thus (33) can be solved with

[TABLE]

satisfying the normalization condition (17).

Next, we assume that the following matricial factorization holds,

[TABLE]

in terms of $2\times 2$ matrices ${\bm{\mathsf{Y}}}_{m}$ to be determined. This factorization is suggested by the results of Coleman and FoxCF63JCP in the fully irreversible regime. Now, substituting the assumption (38) into the equations (34), (II.2),… of the hierarchy, we observe that they can all be solved if the matrices ${\bm{\mathsf{Y}}}_{m}$ satisfy the following relation,

[TABLE]

Introducing the $2\times 2$ matrix

[TABLE]

we see that the solution of Eq. (39) can be written in the form

[TABLE]

Inserting Eq. (41) back into Eq. (40), we obtain the self-consistent matrix equation

[TABLE]

that can be solved by direct numerical iteration for the given values of the rates. Finally, the mean growth velocity (18) is found to be equal to

[TABLE]

with ${\rm tr}$ denoting the trace of $2\times 2$ matrices. In this regard, the matrix (40) is called the velocity matrix. Once this latter is obtained by solving the self-consistent equation (42), we can get the matrices (41), which determine the composition of the copolymer sequences.

Besides, the diffusivity $\cal D$ can be calculated similarly at next order $q^{2}$. According to the central limit theorem, the length probability distribution will thus be given in the long-time limit by the following Gaussian distribution

[TABLE]

Therefore, the solution of the coupled kinetic equations can be expressed for $t\to\infty$ as

[TABLE]

in terms of Eq. (44) and the stationary probability distribution

[TABLE]

to find the sequence $m_{1}\cdots m_{l}$ of length $l$ in the reactive state $i$. The normalization condition

[TABLE]

is satisfied because of Eqs. (17), (31), (II.2),… Moreover, the stationary probability to find the sequence $m_{1}\cdots m_{l}$ of length $l$ in any reactive state is given by the expression

[TABLE]

Since these probabilities cannot be factorized into conditional probabilities as for Markov chains, we conclude that the sequences are nonMarkovian.

All these results can be extended to mechanisms with more than two reactive states, replacing the $2\times 2$ matrices by corresponding $I\times I$ matrices where $I$ is the number of reactive states.

We note that the results of Ref. AG09, are recovered if $w_{i\to j}=0$, as shown in Appendix A.

II.3 Mean, variance, and correlation function

In order to satisfy Eqs. (31), (II.2),…, we have the property that

[TABLE]

implying

[TABLE]

Because of the definition (24), the vector

[TABLE]

is the right-eigenvector of the $2\times 2$ matrix ${\bm{\mathsf{R}}}$ associated with the eigenvalue $\Lambda_{1}=1$. The corresponding left-eigenvector is denoted by the vector $\boldsymbol{\eta}_{1}$. We thus have that

[TABLE]

with a second eigenvalue $|\Lambda_{2}|<1$ and the biorthonormality condition $\boldsymbol{\eta}_{\alpha}^{\rm T}\cdot\boldsymbol{\xi}_{\beta}=\delta_{\alpha\beta}$, where the superscript ${\rm T}$ denotes the transpose.

The mean value of some function $f(m)$ of the monomeric unit $m$ at the location $k$ of the chain is defined by

[TABLE]

Using Eqs. (48) and (49), we get

[TABLE]

where

[TABLE]

so that the mean value is given by

[TABLE]

The variance ${\rm Var}(f)\equiv\langle f^{2}\rangle-\langle f\rangle^{2}$ can be evaluated similarly.

Besides, the correlation function of the function $f$ is defined by

[TABLE]

where

[TABLE]

giving $\langle f^{2}\rangle$ for $n=0$, and

[TABLE]

for $n>0$. Consequently, we have that

[TABLE]

Therefore, the correlation function is given by $\Gamma(0)={\rm Var}(f)$ for $n=0$, and

[TABLE]

for $n>0$. Since $|\Lambda_{2}|<1$, the correlation function decreases exponentially as $\Gamma(n)\sim\exp(-\gamma n)$ for $n\to\infty$ with the rate $\gamma=-\ln|\Lambda_{2}|$. In the following, we shall use the normalized correlation function, $C(n)\equiv\Gamma(n)/\Gamma(0)$.

II.4 Bulk probabilities

The bulk probability $\bar{\mu}(m)$ to find the monomeric unit $m$ anywhere inside the grown copolymer chain can be obtained as the mean value of the indicator function of the monomeric unit $m$ given by the corresponding Kronecker symbol: $f(m_{k})=\sigma_{m}(m_{k})=\delta_{m,m_{k}}$. Therefore, we get

[TABLE]

Similarly, the bulk probability to find the subsequence $m_{1}\cdots m_{k}$ anywhere in the chain is given by

[TABLE]

If we introduce the $2\times 2$ matrix

[TABLE]

Eq. (63) reads

[TABLE]

In general, neither the tip probabilities (48), nor the bulk probabilities (65) have an expression that factorizes as for Bernoulli or Markov chains, which confirms the nonMarkovian character of the macromolecular chains yielded by multistate copolymerization processes.

II.5 Fully irreversible regime

This regime has been studied in Ref. CF63JCP, in the context of polymerization with isotactic and syndiotactic placements, correponding respectively to $m={\rm I}$ and $m={\rm S}$. Fully irreversible regimes are defined by supposing that the detachment rates are negligible in front of the attachment rates:

[TABLE]

Accordingly, the solution of Eq. (42) is given by

[TABLE]

and the mean growth velocity (43) is thus equal to

[TABLE]

which reads

[TABLE]

with the notations $k_{i}=\sum_{m}k_{i,+m}$ ($i=1,2$), $\lambda_{a}=w_{2\to 1}$, $\lambda_{b}=w_{1\to 2}$, and the monomeric concentration $c_{m}=[{\rm M}]$, showing the equivalence with Eqs. (2.4a)-(2.5) of Ref. CF63JCP, . According to Eq. (41), we also have that

[TABLE]

In order to compare with the related results of Ref. CF63JCP, , we introduce the transformation

[TABLE]

such that

[TABLE]

The correspondence with the matrices used in Ref. CF63JCP, is established according to

[TABLE]

for $m\in\{{\rm I},{\rm S}\}$. A direct calculation using MathematicaMathematica shows that

[TABLE]

where

[TABLE]

with $k_{i}=k_{i,+{\rm I}}+k_{i,+{\rm S}}$, which corresponds to Eq. (2.7) of Ref. CF63JCP, , and

[TABLE]

for $m\in\{{\rm I},{\rm S}\}$, which corresponds to Eqs. (3.9a)-(3.9d) of Ref. CF63JCP, . Therefore, the results of Ref. CF63JCP, are precisely recovered in the fully irreversible regime.

III Thermodynamics

III.1 Entropy production

For isothermal-isobaric stochastic processes ruled by kinetic equations

[TABLE]

where $W(\omega\to\omega^{\prime})$ are the transition rates between the coarse-grained states $\omega$ and $\omega^{\prime}$, the link to thermodynamics is established by using the relations

[TABLE]

giving the ratio of the rates of opposite transitions in terms of Gibbs’ free energies $G(\omega)$ associated with the coarse-grained states $\omega$, while $\beta=(k_{\rm B}T)^{-1}$ is the inverse temperature expressed in terms of the temperature $T$ and Boltzmann’s constant $k_{\rm B}$. Gibbs’ free energy $G(\omega)$ is related to the enthalpy $H(\omega)$ and the entropy $S(\omega)$ of the coarse-grained state $\omega$ by $G(\omega)=H(\omega)-TS(\omega)$. Furthermore, the overall thermodynamic entropy of the statistical sample described by the probability distribution $P_{t}(\omega)$ is defined by

[TABLE]

In general, the time derivative $dS_{t}/dt$ of the entropy can be separated into the rate of entropy exchange with the environment of the growing copolymer

[TABLE]

expressed in terms of the mean enthalpy $\langle H\rangle_{t}=\sum_{\omega}P_{t}(\omega)\,H(\omega)$, and the rate of entropy production

[TABLE]

which is always non-negative in accordance with the second law of thermodynamics.

These considerations apply in particular to the kinetic equations (4) where the coarse-grained state $\omega$ is defined by the sequence $m_{1}\cdots m_{l}$ in the reactive state $i$.

In the regime of steady growth where the mean growth velocity is positive $v>0$, the entropy production rate (85) can be evaluated using the probability distribution (45) that is the solution of the kinetic equations in the long-time limit.AG09 ; AG08 On the one hand, the mean value of the Gibbs free energy is given by $\langle G\rangle_{t}\simeq\langle l\rangle_{t}g$ in terms of the mean length $\langle l\rangle_{t}\simeq vt$, so that

[TABLE]

with the mean Gibbs free energy per monomeric unit

[TABLE]

On the other hand, Eq. (45) also implies that

[TABLE]

with the sequence disorder per monomeric unit

[TABLE]

which is always a non-negative quantity, $D\geq 0$. Therefore, the entropy production rate can be expressed as

[TABLE]

with the dimensionless entropy production per monomeric unit, called the affinity,

[TABLE]

the dimensionless free-energy driving force defined in terms of Eq. (87) and the thermal energy $k_{\rm B}T$ as

[TABLE]

and the sequence disorder defined by Eq. (89).AG09 ; AG08

Computing the time derivative of the mean Gibbs free energy $\langle G\rangle_{t}$ with the kinetic equations (81), we get

[TABLE]

where we used Eq. (82). According to Eqs. (45) and (86), we find that the free-energy driving force is here given by

[TABLE]

We notice that the velocity is similarly expressed by Eq. (18).

If the detachment rates are negligible in front of the attachment rates, the free-energy driving force (94) is arbitrarily large, so that the entropy production rate becomes infinite in the fully irreversible regime.

III.2 The equilibrium limit

The state of thermodynamic equilibrium is identified by the principle of detailed balance, which requires that

[TABLE]

This imposes the following constraints on the attachment and detachment rates:

[TABLE]

for all $m=1,2,...,M$. If the attachment rates are proportional to the corresponding monomeric concentrations according to the mass action law, we thus have that

[TABLE]

for all $m=1,2,...,M$, independently of the reactive state $i$. These constraints must be satisfied for the existence of a thermodynamic equilibrium limit.

Now, taking $l=1$ and $m=m_{l}$ in Eqs. (96)-(97), and using Eq. (99), we get

[TABLE]

for $i=1,2$. Therefore, the equilibrium tip probabilities are given by

[TABLE]

because of Eq. (17). Since the probability distribution (101) should also be normalized to unity, we must have

[TABLE]

Moreover, Eq. (38) for $r=1$ gives $\psi(m,i)=\sum_{j}\left({\bm{\mathsf{Y}}}_{m}\right)_{ij}\psi(j)$, while the probabilities $\psi(i)$ define the vector (51), so that

[TABLE]

As a consequence, the sequence probability (48) factorizes as

[TABLE]

showing that the sequences form a Bernoulli chain in the special limit where the equilibrium conditions (95)-(97) are satisfied. Because of Eq. (62), the equilibrium bulk probabilities of the monomeric units are also given by

[TABLE]

Furthermore, the mean growth velocity is equal to zero in the equilibrium limit. Indeed, Eq. (99) has the matrix form ${\bm{\mathsf{W}}}_{+m}={\rm e}^{\varepsilon_{m}}{\bm{\mathsf{W}}}_{-m}$, so that Eq. (40) implies that ${\bm{\mathsf{V}}}\cdot\boldsymbol{\Psi}=0$ at equilibrium and the mean growth velocity (43) is thus vanishing. Accordingly, we find that the equilibrium free-energy driving force (94) can be expressed as

[TABLE]

in terms of the equilibrium bulk probabilities (105). The equilibrium value of the affinity (91) is also equal to zero, which implies that the equilibrium value of the sequence disorder is given by $D_{\rm eq}=-\varepsilon_{\rm eq}$, so that

[TABLE]

in agreement with the result that the copolymer is a Bernoulli chain in the equilibrium limit. As required, the entropy production rate is thus equal to zero at equilibrium.

Besides, the correlation function (61) reduces to $\Gamma(n)={\rm Var}(f)\,\delta_{n,0}$ in the equilibrium limit. Indeed, we have that $\boldsymbol{\eta}_{2}^{\rm T}\cdot{\bm{\mathsf{F}}}\cdot\boldsymbol{\xi}_{1}=\langle f\rangle_{\rm eq}\,\boldsymbol{\eta}_{2}^{\rm T}\cdot\boldsymbol{\xi}_{1}=0$ where $\langle f\rangle_{\rm eq}=\sum_{m}f(m)\mu_{\rm eq}(m)$, because of Eq. (103) and the biorthonormality condition $\boldsymbol{\eta}_{2}^{\rm T}\cdot\boldsymbol{\xi}_{1}=0$ between the left- and right-eigenvectors associated with different eigenvalues. There is thus no statistical correlation in the sequences, as expected since they form Bernoulli chains in the equilibrium limit.

IV Illustrative examples

In the following examples, $M=2$ monomeric species A for $m=1$ and B for $m=2$ are considered. The theoretical predictions are calculated with the methods of Appendix B and the growth is simulated using Gillespie’s algorithm,G76 ; G77 as described in Appendix C.

IV.1 Example with thermodynamic equilibrium

The parameters of this example are taken as

[TABLE]

Here, the concentration $c_{\rm A}$ is the nonequilibrium control parameter. This set of parameter values satisfies the condition (99) for the existence of an equilibrium limit at $c_{\rm A}=1/6$. Accordingly, the thermodynamic equilibrium limit is reached at the concentration where the mean growth velocity vanishes.

For this example, the tip probabilities (46) have been computed using theory and simulations. The results are compared at the concentration value $c_{\rm A}=5$ in Table 1, showing the tip probabilities $\psi(i)$, $\psi(m_{l},i)$, $\psi(m_{l-1}m_{l},i)$, and $\psi(m_{l-2}m_{l-1}m_{l},i)$. On the one hand, these probabilities are obtained in theory by solving Eq. (42) by iterations to get the velocity matrix that is next used to find the matrices (41) and thus the tip probabilities (46). On the other hand, the growth is simulated during the time lapse $t=100$ with Gillespie’s algorithm to generate a sample of $10^{7}$ sequences. With this sample, the probabilities $p_{t}(i,l)$, $p_{t}(m_{l},i,l)$, $p_{t}(m_{l-1}m_{l},i,l)$, and $p_{t}(m_{l-2}m_{l-1}m_{l},i,l)$ are first computed by statistics to next obtain the tip probabilities by summing over the different values of the length $l$ to get $\psi(i)=\sum_{l}p_{t}(i,l)$, $\psi(m_{l},i)=\sum_{l}p_{t}(m_{l},i,l)$, etc… In Table 1, we see the excellent agreement between theory and simulations in this example. The nonMarkovian character of these probability distribution is also confirmed by these results. For the same data, Fig. 1 shows the probabilities $p_{t}(i,l)$ and $p_{t}(m_{l},i,l)$ as a function of the length $l$ in comparison with the theoretical prediction (44)-(45) that the probability distributions should be Gaussian after a long enough time according to the central limit theorem. Here also, agreement is observed between theory (lines) and simulations (data points).

Furthermore, the different quantities of interest have been investigated as a function of the concentration $c_{\rm A}$. Figure 2 depicts the mean growth velocity (43) and the bulk probability (62) to find the monomeric unit A anywhere inside the grown chain versus the concentration $c_{\rm A}$. Again, the data points show the simulation results with Gillespie’s algorithm and the lines the theoretical predictions. As the concentration $c_{\rm A}$ increases, the composition of the copolymer in monomeric units A also increases, as seen in Fig. 2.

The thermodynamic quantities are shown as a function of the concentration $c_{\rm A}$ in Fig. 3. In simulations, the free-energy driving force is computed by adding together the contributions of every jump, as explained in Appendix C. In theory, it is given by Eq. (94). The sequence disorder is obtained with Eqs. (114)-(115). Next, the affinity is calculated by Eq. (91) and the entropy production rate by Eq. (90) in units where $k_{\rm B}=1$.

In this example, the velocity is vanishing at the equilibrium concentration $c_{\rm A}=1/6$. At this special value where the equilibrium detailed balance conditions (95)-(97) are satisfied, the copolymer sequence forms a Bernoulli chain according to Eq. (104) where the tip probabilities (101) take the values $\mu({\rm A})={\rm e}^{\varepsilon_{\rm A}}=1/3$ and $\mu({\rm B})={\rm e}^{\varepsilon_{\rm B}}=2/3$, since $\varepsilon_{\rm A}=\ln(1/3)$ and $\varepsilon_{\rm B}=\ln(2/3)$ for the parameter values (108). As a consequence of Eq. (105), the bulk probabilities take the same values, which is confirmed by the value $\bar{\mu}({\rm A})=1/3$ observed in Fig. 2 at the concentration $c_{\rm A}=1/6$ where $v=0$. Therefore, the equilibrium values of the free-energy driving force (106) and the sequence disorder (107) are here equal to $D_{\rm eq}=-\varepsilon_{\rm eq}=\ln(3/2^{2/3})=0.63651$. The affinity and the entropy production rate are thus vanishing at the concentration $c_{\rm A}=1/6$, as required if equilibrium is reached.

In Fig. 3, the free-energy driving force is vanishing ($\varepsilon=0$) at the concentration $c_{\rm A}=0.70485$, where $v=2.1147$, $A=D=0.68956$, and $d_{\rm i}S/dt=vD=1.4582$. At this concentration, entropy is only produced due to sequence disorder. In the concentration range $1/6<c_{\rm A}<0.70485$, the free-energy driving force is negative, so that the copolymer growth is driven by the entropic effect of sequence disorder. For $c_{\rm A}>0.70485$, the free-energy driving force is positive and the growth is driven by free energy. For large values of the concentration, the free-energy driving force becomes dominant over sequence disorder: $\varepsilon\gg D$. In this strongly irreversible regime where the attachment rates dominate, the entropy production rate can be evaluated with Eq. (94) as

[TABLE]

in units where $k_{\rm B}=1$. For $c_{\rm A}\gg c_{\rm B}$ in the example (108), the entropy production rate increases as $d_{\rm i}S/dt\simeq(8c_{\rm A}/3)\ln(2c_{\rm A})$ with the concentration $c_{\rm A}$, as seen in Fig. 3.

Figure 4 compares the normalized correlation function $C(n)=\Gamma(n)/\Gamma(0)$ of the quantity

[TABLE]

calculated in theory with Eq. (61) in terms of the eigenvalue $\Lambda_{2}$ of the matrix (49) and in simulations with Eq. (120), showing good agreement. For $n=0$, the normalized correlation function is equal to unity. For $n>0$, the correlation function drops to significantly lower values and it decays exponentially at the rate $\gamma=-\ln\Lambda_{2}$. As observed in Fig. 4, the correlation function decays faster and faster as the concentration $c_{\rm A}$ decreases, which is consistent with the absence of statistical correlations in the equilibrium limit where the sequences form Bernoulli chains.

IV.2 Example without thermodynamic equilibrium

The parameters of this example are taken as

[TABLE]

This other set of parameter values is not compatible with the existence of an equilibrium limit.

For this example, the mean growth velocity (43) and the bulk probability (62) for the monomeric unit A are depicted in Fig. 5 as functions of the concentration $c_{\rm A}$ and the thermodynamic quantities in Fig. 6. Again, there is excellent agreement between the simulation results (data points) and the theoretical predictions (lines). Here, the mean growth velocity vanishes at the critical concentration $c_{\rm A}=0.94233$. However, the system remains out of equilibrium, as seen in Fig. 6 showing that, at the critical concentration, the entropy production rate remains positive at the value $d_{\rm i}S/dt\simeq 0.39736$, while both the affinity $A$ and the free-energy driving force $\varepsilon$ diverge as $A\simeq\varepsilon\simeq 0.39736/v$ in consistency with the vanishing of the mean growth velocity. At the critical concentration, the sequence disorder takes the value $D\simeq 0.68$. Since the free-energy driving force remains positive in this example, the growth is always driven by free energy and there is here no regime of disorder-driven growth.

As in the previous example, the free-energy driving force becomes dominant over sequence disorder, $\varepsilon\gg D$, for large values of the concentration $c_{\rm A}$ where the attachment rates dominate. In this strongly irreversible regime, the entropy production rate increases with the concentration $c_{\rm A}$ as $d_{\rm i}S/dt\simeq 2.5\,c_{\rm A}\ln(0.524\,c_{\rm A})$ according to Eq. (109), as observed in Fig. 6.

For the example (111), the normalized correlation function $C(n)=\Gamma(n)/\Gamma(0)$ of the quantity (110) is shown in Fig. 7, where the data points depict the values of Eq. (120) for the simulations and the lines the theoretical functions calculated with Eq. (61). Here also, there is good agreement between theory and simulations with Gillespie’s algorithm. Again, the decay rate is slower at larger than smaller values of the concentration $c_{\rm A}$.

V Conclusion and perspectives

In this paper, theory is developed for the kinetics and thermodynamics of multistate reversible copolymerization processes, in which the growing copolymer may undergo transitions between several reactive states controlling the monomeric attachment and detachment rates and these rates are supposed to be independent of previously incorporated monomeric units. In this case, if the copolymer remained in a single reactive state, its sequences would form a Bernoulli chain. However, as a consequence of the transitions between several reactive states, the growing copolymer sequences instead form a nonMarkovian chain. The probability distribution of the sequences are given in terms of products of matrices associated with every monomeric species and of size equal to the number of reactive states in the mechanism. In general, the chain is nonMarkovian because the matrix products cannot be factorized as for Bernoulli or Markov chains. The matrices of the theory also determine the mean growth velocity, the bulk probabilities of monomeric subsequences, the statistical correlation functions along the sequences, as well as the thermodynamic quantities. In particular, the entropy production rate is shown to be given in terms of the mean growth velocity, the free-energy driving force, and the sequence disorder, confirming the results of Ref. AG08, .

Two illustrative examples are used to compare theory with simulations using Gillespie’s algorithm. For the first example, the rate constants are compatible with the existence of a thermodynamic equilibrium limit. In this example, the entropy production rate vanishes together with the mean growth velocity at the equilibrium value of the tuned monomeric concentration. At this marginal concentration where the equilibrium detailed balance conditions are satisfied, the chain reduces to a Bernoulli chain. Close to equilibrium, there exists a regime where the growth is driven by the entropy effect of sequence disorder, as predicted in Ref. AG08, . For the second example, the rate constants are not compatible with the existence of an equilibrium limit. Accordingly, the entropy production rate keeps a positive value at the critical concentration where the mean growth velocity vanishes and the free-energy driving force diverges together with the affinity. In both examples, the correlation functions are observed to decay exponentially, as predicted by theory. Excellent agreement is found between theory and simulations using Gillespie’s algorithm.

All these results also concern multistate reversible homopolymerization processes yielding sequences with different tacticity. Accordingly, they show that the theory of Coleman and FoxCF63JCP can be extended from fully irreversible to reversible multistate mechanisms, allowing us to include the detachment of monomers beside their attachment in our models. The present theory of multistate reversible copolymerization under low conversion conditions can be further extended to ultimate and penultimate multistate mechanisms where the attachment and detachment rates also depend on previously incorporated monomeric units. For such other kinetics, the sequences form Markov chains instead of Bernoulli chains,GA14 ; G16 if the copolymer remains in a single reactive state, but nonMarkovian chains if transitions occur between several reactive states. The extension of the theory to template-directed copolymerization can also be considered.

Acknowledgements.

This research is financially supported by the Université libre de Bruxelles (ULB) and the Fonds de la Recherche Scientifique -FNRS under the Grant PDR T.0094.16 for the project “SYMSTATPHYS”.

Appendix A The case of kinetics generating Bernoullian chains

In this appendix, we show that we recover previous results for kinetics generating Bernoulli chains if either the transition rates $w_{1\to 2}$ or $w_{2\to 1}$ is vanishing.AG09

If $w_{1\to 2}=0$, the chain stays in the reactive state $i=1$, so that $\psi(1)=1$ and $\psi(2)=0$. Therefore, the mean growth velocity (43) is equal to $v=v_{11}+v_{21}$. However, in Eq. (42), the matrices ${\bm{\mathsf{W}}}_{-m}$ and ${\bm{\mathsf{W}}}_{+m}$ are diagonal and ${\bm{\mathsf{W}}}_{0}$ is upper triangular. Consequently, the iteration (42) converges towards an upper triangular matrix ${\bm{\mathsf{V}}}$, so that $v_{21}=0$. The mean growth velocity is thus given by the solution $v=v_{11}$ of the self-consistent equation

[TABLE]

with $i=1$, as for the growth of a Bernoulli chain staying in this reactive state.

Similarly, if $w_{2\to 1}=0$, the chain stays in the reactive state $i=2$, so that $\psi(1)=0$ and $\psi(2)=1$. In this case, the mean growth velocity (43) is equal to $v=v_{12}+v_{22}$. In Eq. (42), the matrix ${\bm{\mathsf{W}}}_{0}$ is now lower triangular, so that the iteration (42) converges towards a lower triangular matrix ${\bm{\mathsf{V}}}$ and thus $v_{12}=0$. Accordingly, the mean growth velocity is now given by the solution $v=v_{22}$ of the self-consistent equation (112) with $i=2$, as for the growth of a Bernoulli chain staying in this other reactive state.

Furthermore, if $w_{1\to 2}=w_{2\to 1}=0$, the matrix (19) is equal to zero and the matrix equation (42) reduces to

[TABLE]

which is diagonal and decouples into the two previous self-consistent equations (112) for the mean growth velocities of the two types of copolymers, which are thus growing independently of each other.

Appendix B Numerical methods in theory

If there are $I$ reactive states, the matrices are of size $I\times I$. In the reported examples, we take $I=2$. The central equation of the theoretical framework to be solved is Eq. (42) for the velocity matrix. Its solution is obtained by iteration starting from an initial positive matrix ${\bm{\mathsf{V}}}_{0}$. The right-hand side of Eq. (42) gives the next iterate, which is reinserted in the right-hand side, and so on and so forth. This iterative scheme converges towards the solution. In the examples, $10^{4}$ iterations are used to obtain the solution. Thereafter, the mean growth velocity can be calculated with Eq. (43) and the $M$ matrices ${\bm{\mathsf{Y}}}_{m}$ with Eq. (41), giving the tip probabilities according to Eq. (46). In order to get the bulk probabilities and the correlation function, the eigenvalues and the eigenvectors of the matrix (49) are calculated. The bulk probability $\bar{\mu}({\rm A})$ is thus obtained with Eq. (62).

The free-energy driving force is given by Eq. (94), the sequence disorder by Eq. (89), the affinity by Eq. (91), and the entropy production rate by Eq. (90) in units where $k_{\rm B}=1$. In order to obtain the sequence disorder (89), the Shannon sequence entropies

[TABLE]

are first calculated from the tip probabilities (46). The sequence disorder is thus given in principle by

[TABLE]

and, in practice, with the value $D\simeq D^{(4)}-D^{(3)}$.

Appendix C Numerical methods in simulations

The kinetic equations (4) define a Markov jump stochastic process that can be exactly simulated using Gillespie’s algorithm.G76 ; G77 In the examples, we consider the growth of copolymer chains with $I=2$ reactive states and composed of $M=2$ monomeric units. If the current state $\omega$ corresponds to the sequence $m_{1}\cdots m_{l}$ of length $l$ in the reactive state $i$, there are four possible random transitions to a new state $\omega^{\prime}$ that may occur: (1) the attachment of the unit $m_{l+1}={\rm A}$ at the rate $w_{i,+{\rm A}}$; (2) the attachment of the unit $m_{l+1}={\rm B}$ at the rate $w_{i,+{\rm B}}$; (3) the detachment of the unit $m_{l}$ at the rate $w_{i,-m_{l}}$; and (4) the change of the reactive state from $i$ to $j$ at the rate $w_{i\to j}$. The random waiting time $\Delta t(\omega\to\omega^{\prime})$ before the jump is given by an exponential probability distribution of mean value $\tau=\left(w_{i,+{\rm A}}+w_{i,+{\rm B}}+w_{i,-m_{l}}+w_{i\to j}\right)^{-1}$. The transition $\omega\to\omega^{\prime}$ occurs with the probability $P(\omega\to\omega^{\prime})=\tau W(\omega\to\omega^{\prime})$, implying new values for the reactive state, the time, the length, and the cumulated free-energy driving force:

[TABLE]

After many jumps and the elapsed time reaching the value $t$, the mean growth velocity and the free-energy driving force are respectively evaluated by $v\simeq l/t$ and $\varepsilon\simeq E/t$. The determination of the sequence disorder requires the computation of the tip probabilities $\psi(m_{1}\cdots m_{l},i)$ with a large enough statistics, as specified in the figure captions. In the simulations, the tip probabilities are obtained for $l=1,2,3$ and the value of sequence disorder is approximated by $D\simeq D^{(3)}-D^{(2)}$ in terms of the corresponding Shannon sequence entropies (114). The affinity is thus calculated with Eq. (91) and the entropy production rate with Eq. (90). The correlation function is obtained for the function (110) by using

[TABLE]

where $L$ is the total length of the long sequence used in the computation and $\langle f\rangle=(1/L)\sum_{k=1}^{L}f(m_{k})$ is the mean value of the function $f(m)$.