A deterministic model for positional gradients in copolymers

(cid:1) Mathematical model that captures compositional proﬁles of copolymer chains. (cid:1) The model aids in designing polymer chains closest to a speciﬁed composition pattern. (cid:1) Quality measure indicates deviations from the ideal gradient standard. A deterministic modelling approach is developed to predict the internal structure of gradient copolymer chains. A key innovation of the modelling approach is the introduction of a positional variable that gives direct access to quantitative gradient characteristics: the ensemble average composition and the gradient deviation. This positional variable is used to develop multi-dimensional population balance equations that can be solved numerically to calculate gradient quality measures. The methodology is illustrated using the gradient copolymerisation of ethylene and 1-octene via coordinative chain transfer mechanism, which is representative of a variety of polymerisation schemes for gradient copolymers. Simulation results are validated with those obtained by stochastic simulations which, until now, were the only means of predicting detailed gradient quality. (cid:1) 2017 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Living polymerisation systems are attractive because they permit careful control of microstructural properties of polymer chains. An important class of materials that can be produced via living polymerisation is gradient copolymers (Pakula and Matyjaszewski, 1996;Beginn, 2008;Li et al., 2015). Gradient copolymers are similar to block copolymers, but without sharp boundaries between chain blocks. The more gradual gradient structure causes different physical properties in comparison to block copolymers, so that a variety of high-added-value products can be made (Pakula and Matyjaszewski, 1996;Beginn, 2008;Li et al., 2015Li et al., , 2014bZhao et al., 2011;Li and Luo, 2014;Zhou et al., 2013). Gradient copolymers have been produced using different radical polymerisation mechanisms including anionic, cationic, and controlled radical polymerisation methods (Beginn, 2008;Li et al., 2015). For example, the group of Li and Zhu used a modelbased approach with both RAFT (reversible addition fragmentation transfer) and ATRP (atom transfer radical polymerization) to design semi-batch feeding strategies that produce copolymer chains designed to have uniform, linear and hyperbolic gradients (Sun et al., 2007(Sun et al., , 2008Li et al., 2014a;Zhao et al., 2009). Their deterministic models, which were used to select the appropriate reactor operating conditions relied on the method of moments to provide reasonably accurate predictions of M n , polydispersity, and conversion versus time behaviour for the comonomers. Similarly, Ye and Schork developed a method-of-moments model to predict changes in M n ; M w and copolymer composition, as well as average sequence length during semi-batch RAFT and nitroxide mediated radical copolymerisations Schork, 2009, 2010). Their RAFT model was then used to develop optimal semi-batch feeding policies to produce gradient copolymers with linear composition profiles. The models from these two research groups account for chain termination (and in some cases chain transfer) reactions that influence molecular weight distribution (MWD) and disrupt gradient quality. However, their deleterious influence on the quality of the gradient copolymers was largely ignored.
In contrast, Fortunatti et al. (2014) developed a more detailed deterministic model to describe nitroxide mediated radical copolymerisation of styrene and a-methylstyrene in semibatch reactors.
Their model, which accounts for termination and chain transfer, uses probability generating functions to predict the joint distribution of the number of styrene and a-methylstyrene units within the copolymer chains (in addition to M n ; M w and instantaneous copolymer composition). This modelling approach provides helpful information about variability within the population of the copolymer chains produced under different operating conditions. For instance, Fortunatti's model can predict the proportion of low-M n styrene-rich copolymer that forms under different operating conditions. Fortunatti's model, however, does not provide information about how the comonomer units are distributed along the copolymer chains. Zhang et al. (2010) developed a deterministic model to describe semi-batch production of olefin gradient copolymers via Coordinative Chain Transfer Polymerization (CCTP). Zhang's model is similar to several of the reversible-deactivated radical copolymerization models described above, in that it accounts for the influence of chain transfer and thermal termination reactions on M n and M w , but does not account for the influence of these side reactions on gradient copolymer quality (i.e., on the position of the comonomer units along individual chains).
One advantage of the deterministic modelling approaches, when compared with Monte Carlo methods, is that they result in straightforward sets of dynamic balance equations to describe the reacting systems of interest. A second advantage is that the computational effort required to obtain numerical solutions is usually very small, compared with Monte Carlo models, especially if special techniques are used to reduce the number of balance equations that needs to be solved. Some examples of such techniques include the method of moments Schork, 2009, 2010;Zhang et al., 2013), projection methods Iedema, 2013, 2014a;Kryven et al., 2015), and parametrisation with degree distribution (Schamboeck et al., 2017;Kryven, 2017aKryven, ,b,c, 2016. A number of research groups have used KMC models to describe gradient copolymer production via CCTP and a variety of reversible-deactivated radical copolymerisation technologies (Pakula and Matyjaszewski, 1996;Van Steenberge et al., 2012;Valente et al., 2013;Mohammadi et al., 2014;Al-Harthi et al., 2009). The main benefit of the KMC approach is that it permits explicit tracking of a large number of individual copolymer chains, making it relatively easy, compared to deterministic methods, x;y monomer that at the position y of a dormant chain of length x; ending with A at catalyst M d;b x;y monomer that at the position y of a dormant chain of length x; ending with B at catalyst M r;a x;y monomer that at the position y of a living chain of length x; ending with A at catalyst M r;b x;y monomer that at the position y of a living chain of length x; ending with B at catalyst p a ðtÞ total concentration of reacted monomers of type A in the polymer P a ðx; y; tÞ the fraction of polymer chains of length x that contain monomer A at position y P a ðy; tÞ probability that a randomly selected chain has monomer A at position y P x dead chain of length x P d;a x dormant chain of length x; ending with A at catalyst P d;b x dormant chain of length x; ending with B at catalyst P r;a x living chain of length x; ending with A at catalyst P r;b x living chain of length x; ending with B at catalyst t end reaction time x chain length y chain position Fortunatti et al. (2014) andD'hooge et al. (2016) to predict the joint molecular weight and composition distribution of gradient copolymers. Sequence length distributions and other measures of gradient quality are also readily computed (Al-Harthi et al., 2009;Van Steenberge et al., 2012;Wang and Broadbelt, 2009a, 2009b. For example, Al-Harthi et al. developed a KMC model for gradient copolymer production using ATRP (Al-Harthi et al., 2009). They accounted for chain transfer and termination reactions and predicted time-varying cumulative diad and triad fractions in the copolymer in addition to chain length and chemical composition distributions. Wang and Broadbelt (2009a, 2009b developed KMC models to investigate polymerisation of styrene/ methyl methacrylate, focusing on radical activation and deactivation and the resulting variation in average sequence lengths along the copolymer chains. They used their model to select timevarying comonomer and nitroxide feed rates to achieve several targeted properties of gradient polymer chains (i.e., a target M n , a target average sequence length of styrene near the styrene-rich end and a target average sequence length of MMA near the MMArich end). No target was specified regarding the sharpness of the gradient near the middle of the chains. Van Steenberge et al. (2012) developed a KMC model for gradient ATRP copolymerisations involving combinations of n-butyl acrylate, methyl methacrylate and styrene, accounting for termination reactions but not chain transfer. After reviewing the gradient copolymerisation literature, they pointed out that development of gradient copolymers is inhibited by a lack of experimental techniques to quantify gradient quality. One of the benefits of gradient copolymer models is that they can predict detailed information about chain microstructure that cannot yet be measured, but may have an important influence on the physical properties. They can also be used to gain an improved understanding of the influence of undesirable chain-transfer and termination reactions so that improved reactor operating conditions can be selected to produce high-quality gradient copolymer in a relatively short reaction time (Van Steenberge et al., 2015;Toloza Porras et al., 2013;Fierens et al., 2016).
In an effort to summarise information about the location of comonomer units along KMC-simulated copolymer chains, Van Steenberge et al. built on ideas of Beginn (2008) to devise a gradient quality measure hGDi. hGDi provides an indication of the average deviation for an individual polymer chain from a ''theoretically perfect linear gradient copolymer of the same length" (Van Steenberge et al., 2012). Van Steenberge et al. computed values of hGDi from reactor simulations whose semi-batch operating conditions were selected to produce nearly perfect and much-lessperfect gradient copolymers, tracking the evolution of hGDi with conversion (Van Steenberge et al., 2012;D'hooge et al., 2014). For situations where a steep gradient between two blocks is desired, Porras et al. introduced analogous block deviation variable hBDi, which can summarise deviations between simulated copolymer chains and perfect block copolymers with specified target lengths for the two blocks (Toloza Porras et al., 2013). More recently, Lemos et al. (2015) and Devlaminck et al. (2017) applied hGDi measure for process design.
In the current article, an approach utilising multi-dimensional population balances is introduced to permit deterministic predictions of detailed information about location of comonomer units along the chains. This approach is then used to develop a gradient quality measure hGQ i that indicates average deviations from any target arrangement of the comonomer units along the chains. In addition, we show how this approach can be used to provide detailed information about the effects of deleterious chain termination and associated initiation of new chains on the quality of the gradient composition. Since the quality variable can be computed using a deterministic approach, the computations does not require long simulation times and complex data storage that are common to KMC simulations (Wang and Broadbelt, 2011a,b). As a result, we anticipate that the proposed approach will be wellsuited for use in model-based product optimisation, process design and scale-up studies.

Reaction mechanism of CCTP and the positional variable
In the present study we introduce a bivariate deterministic population balance scheme that allows the modeller to quantify the shift of compositional properties along the length of the copolymer chains. We rely on the CCTP scheme of Mohammadi et al. (2014) shown in Table 1 to illustrate the proposed deterministic modelling concept of the positional variable. This ethylene/1-octene system was selected because it has been investigated experimentally in previous modelling studies (Arriola et al., 2006;Mohammadi et al., 2014;Zhang et al., 2010Zhang et al., , 2013. The reaction scheme in Table 1 includes chain initiation, propagation, formation of dormant chains (via reaction 4 and 5) and an undesirable termination reaction, which also leads to initiation of new chains. Pseudo-living polymerisation occurs when reaction 4 and reaction 5 are fast compared to the propagation reactions (reaction 3). Reaction 4 involves the virgin coordinative chain transfer agent CTA (i.e., a metal alkyl like diethyl-zinc, or triethyl-aluminium) and in reaction 5 the CTA moiety is transferred between a CTA-capped dormant chain. If sufficient CTA is added to the system, the majority of the polymer chains are dormant at any time. A narrow distribution of polymer chain lengths results when fast initiation (via reaction steps 1 and 2) occurs at the start of a batch or semibatch reaction so that the polymer chains grow intermittently over the duration of the batch (Valente et al., 2013). Gradient copolymers are produced if there is a shift in the comonomer concentration in the reaction medium over the course of the reactor operation, which leads to a shift in comonomer incorporation along the length of individual chains.
The first step in the model development is to reformulate the traditional ethylene/1-octene CCTP scheme from Table 1 in terms of the positional variable y, as shown in Table 2. That is to say, for each polymerised monomer unit of type A we track several properties: x, the length of the master chain that the unit belongs to; y, the position of the monomer unit along the master chain; whether the master chain is dormant or active (using superscripts Table 1 Reaction scheme for coordinative chain transfer polymerization (Mohammadi et al., 2014).

1.
Active site formation C ! ka C Ã

Initiation
Chain transfer to polymeryl CTA P r;m x1 þ P d;n x2 ! kcs P d;m x1 þ P r;n x2 ; n; m 2 fa; bg 6. Thermal termination P r;m x ! kth Px þ C Ã ; m 2 fa; bg d and r, respectively); and also the type of the monomer unit at the master chain end (using superscripts a and b). The positional variable y is somewhat similar to positional variables introduced by Szymanski et al. (2016) and Gaspard and Andrieux (2014) who counted backward from the living end of a polymer chain to compute sequence length information near the living chain end. In addition to the traditional population balances of living P r;a x ; P r;b x , dormant P d;a x ; P d;b x and dead P x chains, see Table 1, we also formulate population balances for monomers at position y of living x;y , and dead M x;y chains. In the reaction scheme of Table 2, P r;a x refers to a living chain with x monomer units having an A unit at the living end, whereas M r;a x;y refers to an A monomer unit at position y on a growing polymer chain of length x that ends with an A unit. P d;a x;y are A monomer units that are at position y of these chains. Dead polymer chains P x are produced by thermal termination (reaction 6 in Table 1) and bear monomer units denoted by M x;y .
The monomer-focused model formulation is somewhat unusual since it refers separately to various monomer units of the same chain. As shown in Fig. 1, each chain denoted by P x , is associated with various type A monomer units denoted by M x;y . For instance, chain 1 in the figure is referred to as P r;b 10 , as it is a B-terminated (superscript b) living chain (superscript r) of length x ¼ 10. It has three A monomer units referred to as M r;b 10;2 ; M r;b 10;3 and M r;b 10;8 because the overall chain length is 10 and the A units are at positions y ¼ 2; 3, and 8, labelled in order of the growth direction. Note that type B units along the chains do not need to be tracked explicitly because all units that are not of type A must be of type B.
Consider the two propagation reactions shown schematically in Fig. 2 where an A monomer unit is added to living chains of length 10 and 12, respectively. The main propagation reaction equation (at the right) indicates an increase in the length of the polymer chains and the formation of a new incorporated monomer unit. Reaction equations showing simultaneous change in properties of all A species incorporated into copolymer chains are shown at the left. For the chain shown at the top of Fig. 2, the chain length changes from x ¼ 10 to x ¼ 11 and the new monomer unit appears at position y ¼ 11. In the monomeric conversion reactions, y values for the three monomer units remain fixed at 2, 4 and 10, respectively, and the overall chain length increases from 10 to 11. In the reaction at the bottom of Fig. 2, where a new A unit is added to a copolymer chain with a terminal B unit, superscripts denoting the type of terminal monomer change from a to b. Note that in the monomeric conversion equations related to propagation in Table 2 we write '+A' or '+B' on both sides of the equation: the rates of conversion of monomeric units depend on comonomer concentrations, but these conversions do not consume additional comonomer.

Model equations
Consider a semi-batch situation wherein all of reactant B is present initially in the reactor and reactant A is fed continuously over the course of the reaction, producing a gradient polymer that is initially rich in B and becomes richer in A as the chains grow. The feed rate of A can be adjusted to produce any desired gradient. For simplicity, in the model and associated simulations presented here we assume that A is fed to maintain a constant partial pressure in the reactor headspace, which leads to a constant concentration of unreacted A in the liquid solution. Assume further that there is sufficient solvent in the reactor and that the amount of A fed over time is sufficiently low so that any volume changes in the liquid phase over time can be neglected. Also assume that reaction 1 in  Table 2 Reaction scheme for coordinative chain transfer polymerisation using variable y to track positional composition.

Active site formation
x;y þ CTA; m 2 fa; bg 5. Chain transfer to polymeryl CTA P r;m x 1 þ P d;n x 2 ! kcs P d;m x 1 þ P r;n x 2 M r;m x1 ;y 1 þ M d;n x2 ;y 2 ! kcs M d;m x1 ;y 1 þ M r;n x2 ;y 2 ; n; m 2 fa; bg 6. Thermal termination The material balance presented in Table 3 is subject to the following initial conditions at t ¼ 0 : The solution of this system of ODEs is readily obtained by numerical integration. Table 4, the balance equations for the populations of polymer chains, P r;a x ; P r;b x ; P d;a x ; P d;b x and M r;a x;y ; M r;b x;y ; M d;a x;y ; M d;b

As shown in
x;y can be formulated as a system of linear differential-difference equations with time-dependent coefficients. Note that Eqs. (1)-(4) are the classical polymeric population balances corresponding to the reaction scheme shown in Table 1, whereas Eqs. (5)-(8) are the monomeric population balances that obey the mechanism presented in Table 2. The initiation and propagation terms in Table 4 are written involving the Kronecker delta function: dðxÞ ¼ 1 for x ¼ 0 and dðxÞ ¼ 0 for x -0. Here, the delta function ensures that the initiation mechanism produces only chains of length x ¼ 1, whereas the propagation mechanism adds an extra monomer at only the terminal position y ¼ x. Consequently, in Eq. (5) in Table 4, the propagation mechanism results in two terms. The first term increments the chain-length of a chain containing a monomer at position y. The second term, containing dðy À xÞ½P r;b xÀ1 , appends a monomer unit at the end of a chain, at position y ¼ x. Naturally, the domain of values for x and y is limited by the chemistry, so that ½M r;m x;y ¼ ½M d;m x;y ¼ 0, for x; y < 1 or y > x, and ½P r;m x;y ¼ ½P d;m x;y ¼ 0, for x; y < 1 or y > x.
Although the population balance equations given in Table 4 distinguish between living, dormant, and dead chains, the overall polymer product is formed by the mixture of all three. Nevertheless, having concentrations for populations of these types of chains allows us to derive various distributions describing the polymer product at time t. For instance, the univariate length distribution for polymer chains is: where P 1 x¼1 kðx; tÞ is the total concentration of all chains. Similarly, the bivariate distribution for reacted A-monomers that reside at position y of polymer chains with length x is given by the sum: lðx; y; tÞ ¼ ½M r;a x;y ðtÞ þ ½M r;b x;y ðtÞ þ ½M d;a x;y ðtÞ þ ½M d;b x;y ðtÞ þ ½M x;y ðtÞ: The brackets denote the fact that the rate of this reaction is determined by the chain concentration, but that these chains are neither consumed nor produced by this reaction. The mass balance (Table 2) will show this clearly. However, finding the univariate concentration distributions ½P r;a x ; ½P r;b x ; ½P d;a x and ½P r;b x is easily performed by solving the univariate (chain length) balance (system shown in Table 1).   Table 4 Kinetic balance equations for populations of polymer chains. d dt ½P r;a x ¼ k i cacÃðtÞdðx À 1Þ À k ab c b ðtÞ½P r;a x þk ba ca½P r;b xÀ1 þ kaaca ½P r;a xÀ1 À ½P r;a x À Á À k cs0 cctaðtÞ þ kcs k da ðtÞ þ k db ðtÞ ð Þ þ k th ð Þ ½ P r;a x þkcs½P d;a x kraðtÞ þ k rb ðtÞ ð Þ ; (1) x kraðtÞ þ k rb ðtÞ ð Þ ; ( (4) d dt ½M r;a x;y ¼ k i cacÃðtÞdðx À 1; y À 1Þ À k ab c b ðtÞ½M r;a x;y þk ba ca ½M r;b xÀ1;y þ dðy À xÞ½P r;b xÀ1 þkaaca ½M r;a xÀ1;y À ½M r;a x;y þ dðy À xÞ½P r;a xÀ1 À k cs0 cctaðtÞ þ kcs k da ðtÞ þ k db ðtÞ ð Þ þ k th ð Þ ½ M r;a x;y þkcs½M d;a x;y kraðtÞ þ k rb ðtÞ ð Þ ; (5) d dt ½M r;b x;y ¼ k ab c b ðtÞ½M r;a xÀ1;y À k ba ca½M r;b x;y þk bb c b ðtÞ ½M r;b xÀ1;y À ½M r;b x;y À k cs0 cctaðtÞ þ kcs k da ðtÞ þ k db ðtÞ ð Þ þ k th ð Þ ½ M r;b x;y þkcs½M d;b x;y kraðtÞ þ k rb ðtÞ ð Þ ; Note that each chain may contain more than one A-monomer unit, and the sum P 1 x¼1 P x y¼1 lðx; y; tÞ ¼ p a ðtÞ is equal to the total concentration of reacted A monomers incorporated into the polymer.
Using the chain-positional variable as a second dimension, along with chain length as the first dimension, requires us to confront two-dimensional numerical problems. For a maximum chain length of N monomer units (e.g., N ¼ 10 3 ), a one-dimensional population balance is essentially an ODE system composed of N differential equations. The two-dimensional mass balance equations capturing both chain length and the positional information for each monomer consist of N 2 of such differential equations. Rather than perform a traditional numerical integration of such a large number of differential equations we utilise an approximation method. The projection onto radial basis functions (Kryven and Iedema, 2013) replaces the whole ODE system with a much smaller one that can nevertheless accurately describe the solution. This method has enabled detailed studies of various polymerisation and coagulation systems within chemical Iedema, 2013, 2014a,b;Kryven et al., 2014) and biological (Kryven et al., 2015) contexts.

Gradient quantification
Beginn noted that phenomenon of gradient composition can be well quantified only when the whole ensemble of copolymer molecules is taken into the account. In his review on gradient copolymers (Beginn, 2008), the concept of the ensemble average of the local chain composition is introduced as follows. Let there be N chains in the ensemble. Indicator function F A;i ðyÞ 2 f0; 1g is set equal to 1 if there is monomer A at position y on the i th chain, and set at zero if monomer B occupies that position. Then is the ensemble average fraction of monomer units of type A at position y. Alternatively, assume one draws a monomer unit uniformly at random from all monomers having position y, then F a ðyÞ h i2½0; 1 can be viewed as the probability that the selected monomer is the A monomer. Van Steenberge et al. (2012) used the cumulative properties of F A;i ðyÞ to define a measure of gradient quality as a distance form a target cumulative distribution. The ensemble average of this deviation, hGDi, is then taken as the measure that expresses to which degree the ensemble features the gradient property. Although F a ðyÞ h iand hGDi have been generated using MC models, until now no gradient measures have been explained or calculated using a deterministic approach. We will show now that the population balance model from the previous section produces enough information to quantify the chain gradient in a purely deterministic fashion. The fraction of polymer chains of length x that contain monomer A at position y is P a ðx; y; tÞ ¼ lðx; y; tÞ kðx; tÞ : From probabilistic perspective, P a ðx; y; tÞ is the conditional probability that a chain contains a monomer A at position y given that chain length is x. Furthermore, one may relax the conditioning on chain length and write the probability that a randomly selected chain has monomer A at position y, P a ðy; tÞ ¼ P 1 x¼1 lðx; y; tÞ P 1 n¼y kðn; tÞ : This ratio has the concentration of reacted monomers of Type A at position y in the numerator and the number of all polymer chains of length y or greater in the denominator. The complementary probability of finding a B-unit at position y is P B;t ðyÞ ¼ 1 À P A;t ðyÞ. These probabilities are the ensemble average compositions introduced by Beginn (2008), P a ðy; tÞ ¼ F a ðyÞ and P b ðy; tÞ ¼ F b ðyÞ, and can be readily computed from our model. Similar to the linear gradient quality, a variable of Van Steenberge et al. (2012), which reflects the deviation from an ideal linear gradient, the quality of the gradient composition can also be estimated from the deterministically calculated distribution lðx; y; tÞ. The ideal linear gradient from B to A in a chain of length x is defined by: P ideal a ðx; yÞ ¼ y x when counted in the growth direction, or P ideal a ðx; x À y þ 1Þ when counted in the contra-growth direction. Defining the ideal gradient concept allows us to calculate the average mean square deviation from the ideal gradient, Here, we calculate the quality for both gradient directions and accept the smaller value (a better value). Note, that under current polymerisation mechanism, the predominant majority of chains will comply with the first choice of the direction, that is probability of B-monomer inclusion declines in the direction of chain growth. That said, both directions are implemented in Eq. (11) so as to ensure applicability of hGQ i to polymerisation mechanism in the most general context. When the actual gradient is equal to the ideal one the gradient quality becomes zero, whereas the normalisation factor 3 y ensures that hGQ i is equal to one in the worst-case scenario. In this way, the quality hGQ i is a number between 0 and 1: the ideal gradient quality (linear gradient) corresponds to hGQ i ¼ 0, and pure A or pure B chains correspond to hGQ i ¼ 1. Measure hGQ i can be generalised to an arbitrary nonlinear target profile by replacing the target function and adjusting the normalisation factor.
To summarise, measure hGQ i provides a deterministic alternative to stochastic measure hGDi, although these two measures are not strictly equivalent: hGQ i compares collective profile of the whole ensemble against the target, whereas hGDi is an average of individual deviations of all chains from the target. This difference becomes apparent if one considers an ensemble with large variance of chain length. For example: a mixture of short gradient chains and long gradient chains will score good quality according to hGDi but not according to hGQ i, since the collective compositional profile of many gradient chains with considerably different length might not necessary be gradient. Depending on particular application one may prefer to use one measure over the other.

Numerical experiments and discussion
In this section we present numerical results obtained for the above-introduced population balance model. The reactivity constants were adopted from reference Mohammadi et al. (2014) and are given in Table 5 for reference. Here, we denote ethylene as monomer A and 1-octene as monomer B. The population balance (1)-(8) not only allows computing chain-length distribution, but also to understand how individual monomers are distributed along a polymer chain. From numerical experiments, we are able to deduce the optimal values of input parameters that yield, for target M tar w ¼ 140 kg mol , conversion v tar ¼ 0:93 and c a þ c b;init ¼ 6:7 Á 10 3 mol m 3 , the best gradient quality in the sense of measure (11). In other words, our aim is to solve a constrained minimisation problem, in which the following four input parameters that affect objective function hGQ i are subject to optimisation: (i) initial octene monomer concentrations, c b;init ; (ii) constant ethylene monomer concentration, c a ; (iii) initial concentration of CTA, c cta;init ; (iv) initial concentration of the catalyst, c cat;init .
In order to understand what the parameter space looks like in view of the optimisation constraints, we will discuss them by gradually adding the constraints one-by-one. Firstly, c a þ c b;init ¼ const, and therefore, if the initial total monomer concentration is given, one may replace these two parameters simply by initial monomer fraction, f a ¼ ca ca þc b;init . Secondly, the range of values for final conversion v is not restricted by any other parameter. For this reason, one may fix the initial monomer fraction f a and readily analyse the model for v ¼ v tar . Fig. 3 illustrates how the average molecular weight depends on the values of c cta;init ; c cat;init under these conditions. Such an analysis reveals that M w behaves as a linear function with respect to c cta;init and c cat;init . So that those values of c cta;init and c cat;init that all lead to equivalent, in terms of molecular weight, polymer M w ¼ M tar w , can be parametrised by a straight line, as depicted by the read line in Fig. 3. This result may seem counterintuitive since it implies that M w ¼ M tar w even for c cat;init ¼ 0. However, one must realise that in the limit c cat;init ! 0, the reaction time required to reach the target molecular weight approaches infinity, as depicted in Fig. 4. The behaviour shown in Fig. 3 is expected because the total concentration of living polymer chains in the system is equal to the sum of the initial concentrations of catalyst and CTA. Adding more of either tends to reduce the chain length.
A natural question now arises: do these multiple options of parameters c cta;init þ c cat;init ¼ const, yielding chains with the same average molecular weight, also yield equivalent chains in terms of positional composition and gradient quality? It turns out that these parameter sets all lead to almost equivalent positional profiles having hGQ i ¼ 0:0324 AE 0:0002. So that as long as c cta;init þ c cat;init ¼ const the value of c cta;init separately does not influence the quality of the gradient, but it is still the decisive factor for the reaction time. Taking these considerations into account, we fix c cta;init ¼ 5:767 mol m 3 ; c cat;init ¼ 1:43 Â 10 À3 mol m 3 , as depicted by the withe star in Fig. 3, that leads to reasonable reaction time of t ¼ 1000 s at the target conversion. So that solving the optimisation problem, which expects solving (1)-(8) followed by evaluation of (11) on every iteration of the optimisation algorithm, leads to the optimal-value sets of the input parameters as given in Table 6.   Table 6 Optimal polymerisation parameters as obtained for two distinct sets of constraints.

Constraints
Optimal parameters We will now turn to a more detailed analysis and discuss the compositional profiles of the chains produced in the hGQ ioptimal mode. The optimal gradient profiles P a ðy; tÞ; P b ðy; tÞ, are shown in Fig. 5. In this figure the data obtained from the deterministic population balance model (1)-(8) is validated against the MC simulated data. The corresponding Monte Carlo model was developed in a simple fashion by applying Gillespie's traditional KMC method (Gillespie, 1977). To perform the MC simulation, a small initial reaction volume was specified and the corresponding initial numbers of each type of reactant molecules (i.e., ethylene, 1octene, catalyst and CTA) were computed. Uniformly distributed random numbers were generated repeatedly and used to determine: the duration of the short random time intervals between reactions, they type of reaction that occurred after each interval, and which polymeric molecules were involved in the randomly selected reaction. After each time interval, the number of small molecules of each type was updated as were the chain lengths, compositions and types of the polymer molecules that participated in the selected reaction. This updated composition information was then used to determine probabilities of the different types of reactions that might happen at the end of the next random time step.
In both the KMC and deterministic simulation results in Fig. 5, one notes that P b ðy; tÞ decreases with increasing position from the beginning of the chain (counted in the direction of polymerisation). This feature indeed points at the presence of a statistical gradient in the chain composition. However, this is only part of a larger picture, since P a ðy; tÞ contains information that is obtained by averaging over chains of different length, whereas classes of chains of equal length, may still feature considerable variation in composition profiles.
Exploring eventual differences between chains of one length class has is possible with the two-dimensional population balance model -it allows us to zoom in even further. As shown in Fig. 6a and b, a rapid decay in probability of B inclusion as a function of position can still be observed. The composition patterns, however, clearly depend on the length of a chain. In fact, the chains with the best gradient quality are also the chains close to ensemble average in length. This is not a coincidence but an outcome of the optimisation routine, since the optimisation was carried out for a specific target molecular weight. We also see that the short chains have almost a uniform, B-rich, composition. This observation indicates that creating a narrow chain-length distribution in this polymerisation system is important when high gradient quality has to be achieved. That said, it should be kept in mind that Fig. 6 presents  only a conditional probability: it represents the composition profile of chains of a given length x. The actual chain-length distribution, as given in Fig. 7, reveals that, indeed, the length of the chains is narrowly dispersed around the mean value. This leads to rather a good gradient quality, which besides being expressed by the low value of hGQ i, may also be qualitatively assessed in Fig. 8 that is generated by Monte Carlo simulations. The next question we address is the justifiability of the target conversion. It turns out that increase of the conversion does not improve the gradient quality. As illustrated in Fig. 9 the gradient quality, hGQ i, has a minimum at v ¼ 0:93, which was also the target conversion of the optimisation procedure. Systems that are even closer to full conversation are not necessarily associated with a better quality of the gradient due to excessively long A sequences that appear in the end of the chains. Yet, if the optimisation of the input parameters is reperformed for higher target conversions, the optimal gradient quality improves up to hGQ i ¼ 0:011 at v tar ¼ 0:97, see Table 7.
Finally the last observation regards the effect of the thermal termination mechanism on the gradient quality hGQ i. As demonstrated in Table 8 increased values of termination rate k th are associated to degeneration of the optimal gradient quality. Such a degeneration of hGQ i is caused mainly due to an increase in the variance of the chain-length distribution.

Conclusions
We have presented a mathematical model capable of explaining nuances of gradient copolymerisation in a deterministic way. The central idea of the model is the population balance approach that describes an ensemble of copolymer chains with distributions. Importantly, the model reveals not only how chains with different lengths are distributed but also how individual monomers are distributed along a polymer chain. An in-depth view of the polymer chains has been achieved by introducing the positional dimension: polymer chains are governed by a certain distribution as members of a large ensemble, but these chains, in turn, consist of monomers distributed according to certain statistical pattern along the length of the chains. Until now, such positional information has only been generated by Monte Carlo simulations. The model is the first deterministic model that provides an understanding of compositional structure of individual polymer chains.
The proposed approach could be important when one aims to design copolymer chains with a specific composition pattern. To demonstrate the utility of this chain-positional approach we selected a linear gradient specification rule as the target composition profile. We then optimised the model parameters (i.e. the concentrations of reactants in the semibatch reactor) in order to reach the smallest mean-square deviation from this standard in the final polymer product. The proposed methodology could readily be extended to account for more complicated reactor feed strategies and to optimise reactor operating conditions for other target gradient specifications. Furthermore, the methodology is quite general in that it can be used for other gradient polymerisation mechanisms such as ATRP and NMP. Accounting for termination by combination as an undesirable side reaction is an added complication that will be explored in a future article. Composition profiles of polymer chains play decisive roles in the physical and thermodynamic behaviour of copolymers, e.g. self assembly into regular structures (micelles, brushes) or crystallisation. Consequently, the freedom to design and achieve compositional patterns of polymer chains is essential for fine tuning of such processes. Fig. 9. Gradient quality hGQi, the mean square deviation of the chain composition profile from the linear gradient, versus reaction progress. Gradient quality reaches the best (lowest) value at v ¼ 0:93 but becomes worse as the reaction proceeds further.