Dynamic Networks that Drive the Process of Irreversible Step-Growth Polymerization

Many research fields, reaching from social networks and epidemiology to biology and physics, have experienced great advance from recent developments in random graphs and network theory. In this paper we propose a generic model of step-growth polymerisation as a promising application of the percolation on a directed random graph. This polymerisation process is used to manufacture a broad range of polymeric materials, including: polyesters, polyurethanes, polyamides, and many others. We link features of step-growth polymerisation to the properties of the directed configuration model. In this way, we obtain new analytical expressions describing the polymeric microstructure and compare them to data from experiments and computer simulations. The molecular weight distribution is related to the sizes of connected components, gelation to the emergence of the giant component, and the molecular gyration radii to the Wiener index of these components. A model on this level of generality is instrumental in accelerating the design of new materials and optimizing their properties, as well as it provides a vital link between network science and experimentally observable physics of polymers.


Results
A polymer is a large molecule that consists of many repeat units, the monomers, and is formed as a result of chemical reactions that lead to covalent bonding between the monomers.
Step-growth polymerization does not require an initiator and occurs between monomers that carry reactive functional groups. Many polymers with real world applications are formed as a result of step-growth polymerisation. Figure 1 features a few important examples related to polyesters, polyamides, and polyurethanes. The maximum number of chemical bonds that a single monomer bears is limited by the number and type of functional groups that are present in this monomer. If a system consists of solely two-functional monomers, only linear polymers are formed in the course of the step-growth polymerisation. However, if some (or all) monomers have more than two functional groups, it is possible to form hyperbranched polymers and networks.
In many chemical systems, two monomers bind through an asymmetric reaction that occurs between functional groups of different kinds. When two functional groups of different kinds are reacting, for example, as in the reaction between an acid and an alcohol leading to an ester, we refer to one group as the A-group, and the otheras the B-group. This asymmetric reaction is at the main focus of the current paper. Symmetric reactions occurring between two groups of the same kind, e.g. two alcohols reacting to form an ether, have been covered elsewhere 22 . Figure 1 exemplifies this notation on a few cases of polymers that feature linear and branched topologies. For each case, we indicate the structural formulas of the relevant monomers, highlight the functional groups, and give the corresponding AB notation. In Figure 2 the asymmetric reaction between A and B groups is illustrated on the example of an A 2 monomer (a monomer with two A-groups) reacting with an AB 2 monomer (a monomer with one A-group and two B-groups).
Before introducing the random graph model, we briefly summarise the terminology commonly used in graph theory. A directed graph consists of nodes and directed edges connecting them. A subgraph of a graph, in which any two nodes are connected by an undirected path is called a weakly connected component. In this work, we drop the prefix weakly, and refer to these components as connected components. When representing a polymer system as a graph, the monomers are identified with nodes, the chemical bonds with edges and thus a polymer molecule with a whole connected component. As the two sides of a chemical bond in the AB-reaction are not identical, we represent this asymmetry with directed edges. Without loss of generality, the directionality is defined as pointing from the B-group towards the A-group. The graph representation and the mapping of a reacted A-/B-group to an in-/out-edge is depicted in Figure 2.
The number of bonds that a monomer is bearing equals to the degree of the corresponding node. Generally speaking, monomers with distinct numbers of bonds may have different concentrations. To capture these differences, we refer to the node degree distribution, u(i, j), which defines the probability of a randomly chosen node to have i adjacent in-edges, and j adjacent out-edges, and therefore, u(i, j) is proportional to the concentration monomers with i and j reacted A-and B-groups. Figure 3b, demonstrates that the directed degree distribution u(i, j) and the projected undirected degree distribution = ∑ + = u k ui j ( ) ( , ) i j k that ignores the direction of the edges may lead to different sizes for connected components. In this extreme example, the directed degree distribution only allows connected components of size s = 3, whereas the undirected degree distribution does not limit these sizes at all. It turns out that much of the global information about the polymer system can be deduced from the degree distribution. Hyperbranched molecules appear in a broad range of topologies that result from various chains of reactions, which occur between selected functional groups. One important feature here is the asymmetry of such reactions: they often occur between functional groups of complementary types (as opposed to symmetric reactions occurring between groups of identical types). To capture this variability we employ the configuration model that is defined by a directed degree distribution. The degree distribution reflects the state of the polymer system as driven by the chemical kinetics, and therefore, this distribution is time-dependent. The configuration model maximises entropy of all possible configurations that satisfy a given degree distribution at a time point of interest. This  means that the model is egalitarian with respect to functional groups: every pair of complementary functional groups has equal probability to establish a bond. This principle is illustrated in Figure 3a.
The output of the configuration model is then processed to obtain various global properties of the polymer network, as for instance, the molecular weight distribution (probability that a randomly chosen node belongs to a weakly connected component of size s), the gel fraction (probability that a randomly chosen node belongs to the giant component), the mean-square gyration radii (related to the Wiener index of connected components), and the length of the average shortest path. The mathematical derivations of these results are presented in the Methods.
Reaction kinetics and degree distribution. The evolution of the degree distribution is governed by the reaction kinetics of the step-growth polymerization. This process converts an arbitrary pair of A-and B-groups into a bond, which is represented as an edge between nodes in the model. Since a monomer may carry multiple A-and B-groups, probability that a pair of monomers react is dependent on the number of unreacted functional groups they carry. The formal mechanism for the reaction between two monomers is given by: where vector (i, j, I, J) denotes the state of a monomer: I, J ≥ 0 are the numbers of respectively A-and B-groups on this monomer; i = 0, 1, …, I and j = 0, 1, … J denote the number of groups of respectively type A or B that have already been converted into bonds by the reaction. For each monomer, the indices (I, J) are defined a priori, whereas (i, j) change over time. The reaction rate is given by the product of the rate constant k AB , the number of unreacted A-groups on the first monomer (I − i), and the number of unreacted B-groups on the second monomer (J′ − j′). The following assumptions are made: (1) the reactivity for any pair of A-and B-groups is equal; (2) the reactivity does not change throughout the process. Let M i,j,I,J (t) be the probability that a randomly chosen monomer has configuration (i, j, I, J) at time t. This probability is proportional to the concentration of the monomers. In (a) the colour of the bond indicates which monomer provided the A-and the B-group. In (b), the graph representation, the type of the functional group is stored in the directionality of the edge. An in-edge corresponds to a reacted A-group, an out-edge to a reacted B-group. Nodes are signed "half-edges" according to the bivariate degree distribution that provides an input for the model; (right:) these nodes are connected randomly so that the degree distribution is strictly satisfied. (b) An example of a simple degree distribution and the corresponding ensemble of connected components for the case of undirected and directed graphs. Note that in this example there is a qualitative difference between the structure of undirected and directed chains. The time variation as governed by the process eq. (1) is described by the corresponding master equation, see Methods. In the general case of more than two distinct functional groups (e.g. A-, B-, C-, D-groups) and several reaction mechanisms (e.g. AB-, CD-reaction), the master equation becomes more complex due to a combinatorial number of monomer species of defined type and state. In that case, automated reaction networks can be applied to algorithmically construct the corresponding master equation 23 .
The temporal degree distribution u(i, j, t) is directly deduced from M i,j,I,J (t) by summating over all monomer types I, J. As discussed in Methods, the step-growth process eq. (1) leads to the following degree distribution u(i, j, t): This expression is given by the product of binomial distributions for the in-and out-edges, with p A (t) and p B (t) being the probability that a random A-/B-group is reacted (also referred to as A-/B-conversion), see eqs (12) and (13). The probability distribution P(I, J) defines the initial concentration of monomer types and is referred to as the monomer functionality distribution. Probabilities P(I, J) provide the sole input to the model. In some cases, the degree distribution u(i, j, t) can be measured directly by nuclear magnetic resonance (NMR). Due to the chemical shift in the NMR spectrum, every monomer state has its own distinct frequency. However, the bigger the monomer is and the more functional groups it has, the harder it is to identify the states. In Figure 4 we compare the degree distributions predicted by eq. (2) with the experimental NMR data from ref. 50 , and also, for a different polymerisation system, against MD simulation data from ref. 40 .
Gelation point and gel fraction. The gelation point marks the transition of the system from a liquid-like to a solid-like state during the polymerization process. After this transition, an increasing positive fraction of monomers becomes part of the single gel molecule. The transition is typically observed by measuring the fraction of the insoluble part of the polymer, or performing rheological experiments. From the network theory perspective, this is a well-known phenomenon as the gel transition point corresponds to the emergence of the giant component in the network topology. The size of the giant component is of the order of the whole system size, and we therefore quantify the gel size g f as the probability that a randomly chosen monomer belongs to the giant component. The Methods section explains how g f can be calculated if the degree distribution is known.
The above described process of polymerisation is closely related to percolation on networks. The phenomenon of percolation is always studied on a specific networks, e.g square lattice, Bethe lattice, or any arbitrary network. In polymer chemistry, it is mostly introduced as a process of randomly adding bonds to an empty template of the network 51 . This process is precisely reverse to the process of removing bonds from the same full network. Let p be the probability that a randomly selected edge is not removed.
Under this notation, percolation can be thought of as a temporal process that starts at p = 0 and ends at p = 1. Moreover, this process is known to feature a phase transition at critical probability p critical , the point at which the giant component appears. When there are equal amounts of A and B functional groups, this process is precisely reverse to the step-growth polymerisation where the edges are being added to a network randomly, and p critical coincides with the gel-point conversion of functional groups. If there are more B-groups than A-groups present in the system initially, the A-groups are limiting the reaction and we conventionally refer to the conversion of A groups as the conversion: p = p A . Without loss of generality it may be assumed that the number of A-groups is less or equal to the number of B-groups. The moment in time when gelation occurs is completely defined by the proportion of monomers of different functionalities that are present initially in the system. Generally speaking, the higher the functionality, the earlier gelation occurs in time and/or conversion, and a precise quantitative estimate of the gelation conversion is discussed in Methods. It turns out that one can connect the critical conversion directly to the functionality distribution P(I, J): . It is important to note that in the case of a symmetric functionality distribution P(I, J) = P(J, I) (see also Methods), eq. (3) gives the same result as its counterpart for an AA-system (undirected system 22 In the case of P(K) being distinct from zero for only one K-functional species, eq. (3) simplifies to the Flory-Stockmayer equation for gelation = − p z critical 1 1 , with z = I + J the number of first neighbours. See Table 1 for illustration.
Equation (3) allows one to screen a vast number of systems and determine their gel-point conversions if such occur. We also can deduce from the theory discussed in Methods that some systems will never form gel. Here again, one can identify whether a system of given monomer functionalities gelates by studying P(I, J). Namely, a monomer system forms gel if the following inequality is satisfied: For the physical properties of the final material, both factors play a definitive role: when does gelation start and how does the gel fraction evolve in the course of the polymerisation? The growth rate of the gel fraction, g f , is determined by the functionalities of the initial monomers and their concentration distribution. In Figure 5, a few examples illustrate different types of behaviour of the gel buildup. In these examples, we optimise the initial functionalities and concentrations of monomers to reach two final target properties: (1) a fixed gel point conversion of either p A,critical = 0.5, as depicted by the solid lines, or p A,critical = 0.33, as indicated by the dashed lines; (2) we distinguish three different types of growth behaviour: (a) a steep growth with most monomers being incorporated into the gel rapidly after gel point, (b) a slow growth with the gel reaching full size only at full conversion, (c) a slow growth with the gel never reaching the system size. Behaviour (a) is observed for systems with purely high-functional monomers, (b) for systems with few high-functional monomers and many 2-functional monomers, and (c) for few high-functional monomers and many 1-functional monomers that act as terminal units. The reason for the gel in (c) never reaching the full system size is the formation of small connected components that stop growing because of having all functional groups being capped with one-functional terminal units. For example, when a component is composed of one 6-functional monomer connected to six 1-functional monomers. The Methods section gives the general equations for the gel-point conversion and gel fractions.  These results are also interesting when studying polymer ageing and degradation. Consider a degradation process under which every chemical bond dissociates independently with equal probability. This process is reverse to the introduced polymerization process, and the gel fraction is a measure of how strongly the system is interconnected. Clearly, the systems of type (a) will show a different behaviour during degradation than systems of type (b). For (a), the system will stay connected for a long time, but will eventually collapse into many small pieces quite abruptly. Type (b) systems will show a more continuous, and therefore more predictable, degradation behaviour, which might be more desirable from application point of view. Molecular size distribution, averages and asymptotes. From a network theory perspective, a separate polymer molecule is a connected component. The sizes of the latter are typically characterised by a size distribution. There are two common ways that such distributions can be defined. The molecular weight distribution w(s) corresponds to the probability that a randomly chosen monomer belongs to a connected component of size s, whereas the molecular size distribution n(s) is the probability that a randomly chosen component has size s. One can be converted into the other by an appropriate weighting and normalisation: n(s) = Cs −1 w(s). In Methods we present an exact equation that connects the degree distribution with the molecular weight distribution w(s). As a general rule, the exact values of w(s) can be computed spending  s s ( log ) 2 multiplicative operations, and in a special case of only one initial monomer type, the analytic expression for the molecular weight distribution is given in Methods.

Directed Undirected Flory-Stockmayer
The global behaviour that is observed in all polymerising systems can be summarised as follows: Initially, all monomers are unconnected, thus only molecules of size s = 1 are present. With increasing conversion p A , larger molecules emerge. The size distribution features the exponential decrease at the tail, and becomes broader with progressing conversion until the gel point p A,critical is reached. Only at this single point the size distribution becomes scale-free. Figure 6a demonstrates this behaviour on an example. After the gel transition point, the size distribution describes only the soluble part of the system, so that the size of the gel is given by the gel fraction . Furthermore, the size distribution returns to its exponential behaviour and becomes narrower with increasing conversion.
Surprisingly, there are two distinct types of polymerisation systems featuring different types of asymptotic behaviour. Most of polymer systems feature a size distribution with asymptote This, for instance, includes the A 2 :B 3 = 3:2 system as illustrated in Figure 6a. However, some polymers may also feature a different asymptotic mode, namely This asymptote arises in all systems of type AB n , for n > 1. These two distinct types of asymptotic behaviours can also be attributed to the different types of kernel functions (additive versus product kernel) in the underlying aggregation equation 11,12 . In Figure 6b, this peculiar case of asymptotic behaviour is illustrated by comparing two very similar polymer systems AB 2 and A 2 + B 3 that yet feature different asymptotic modes. The gel transition is also noticeable in the evolution of weight average molecular weight M w , which features a singularity at the critical point. The evolution of M w as predicted by the theory is compared against stochastic simulations in Figure 6c. The figure shows good agreement between the theory and the numerical simulation except for critical conversion. At this point, the stochastic simulations (scatter plot) suffer form the small-system-size effect. The Methods section gives analytical equations for the weight average molecular weight in the pre-gel and gel regimes.
Interestingly, in some cases molecular size distributions feature oscillations. One such example is given in Figure 6d, depicting the system AB 2 :A 1 :B 1 = 1:2:1. At full conversion p A = 1 (dark red line) only molecules of specific favoured sizes are present. At lower conversions, also other sizes occur that exhibit strongly reduced probability as compared to the favoured sizes. In this example the concentration of odd-size component is favoured, however, depending on the distribution of functionalities such oscillations may occur with arbitrary large "periods". It turns out, that the monomers of functionality one play an important role in these oscillations, as they terminate the growth of polymer molecules and thus fix sizes of these molecules at a constant value.
Gyration radius. Consider a branched polymer molecule that is composed of s monomers. The actual volume this molecule spans is related to how'branched' this molecule is. In systems that contain no gel, or are below the gel transition, it is conventional to characterise this volume by the a quantity called gyration radius R g (s), which can also be estimated by light scattering experiments in a polymer solution 52 .
Linear chains feature , where b is the Kuhn length. In Methods, we derive the analytical equation that links the degree distribution and the mean square gyration radius for  s 1. Figure 7a shows how one can influence R g (s) by tuning the set of initial monomers. This figure also compares the theoretical gyration radii against gyration radii obtained form stochastically generated networks. An alternative way of looking at the gyration radius is the contraction factor 53 , which is defined by . The contraction factor tells us how much more compact the actual molecule is in comparison to the linear one having the same number of monomers. Figure 7b compares the theoretical predictions of this quantity versus simulations.
Another unexpected result that is revealed by directed random graph theory, is that during the progress of step-growth polymerisation, for all molecules of size s the mean gyration radius is constant in time. Since the size distribution does change in time, this time-independency is lost if one calculates an average of this quantity over different molecular sizes. It is likely that other than step-growth polymerisation processes do not feature time-invariant gyration radii.
The gyration radius is also directly proportional to the Wiener index 54   Scaling of the node neighbourhood size and criticism of well-mixing assumption. Since the radius of gyration characterises only finite-sized molecules rather than the gel, which is virtually infinite in size, we are in need of devising an additional measure for the degree of connectedness of the gel. With this aim, we investigate the number of nodes at distance l from a randomly selected one, which we shortly refer to as  l . Here, we utilise the topological notion of distance: nodes i and j are at distance l if the shortest undirected path connecting these nodes has length l. Since the giant component is infinite in size, the larger the distance l the more nodes are incorporated in the volume of a sphere. In fact, we are mainly interested in the way this quantity asymptotically depends on  l 1. Newman derived an expression for a similar quantity for directed paths 41 , whereas in the polymer context we are interested in the weak sense of connectivity. This means that we consider the shortest paths that does not respect the directionality of the bonds.
An analytical expression for this behaviour is derived in Methods. We observe that the average path length in the gel exponentially depends on l: l l l / 0 where C and l 0 are constants. Note that structures that reside in three-dimensional Euclidean space must feature a different scaling law: where 1 ≤ d f ≤ 3 is the cluster growth dimension. The estimate given in eq. (7) is quite a discouraging result, as a network that features an exponential scaling cannot be physically embedded in the Euclidean space under the condition that the nodes are uniformly distributed with constant density. This unphysical exponential growth of a node neighbourhood is caused by one of the fundamental assumptions made in the random graph model, and many other popular models that do not explicitly track the spatial configuration of the network also suffer from the same criticism. In fact, we refer here to one of the most commonly used assumptions of a well-mixed system: any two functional groups react with an equal probability irrespectively of their location in the topology. In real systems, however, monomers interact with the rest of the network that may locally hinder them to react. That said, it is important to note, that the scaling given by eq. (7) does not hold before, and precisely at, the gel transition, and the well-mixed system assumption may remain a good approximation in these regimes.

Discussion
This paper employs recent developments in network science to formalise the step-growth polymerisation process as a problem fully described by network generation. In order to do so, we proposed to view the polymer architecture as a directed network that is described by a dynamic degree distribution. Although the physical connection between monomers by means of covalent bonding is completely symmetrical, the directionality of the edges keeps record of the asymmetry of the chemical reaction that created the corresponding covalent bond. This approach allows a classification and a general treatment of a vast range of real-world polymerisation problems and can be used for optimisation and design of new materials. As a general rule, the parameters of step-growth polymerising systems comprise a high-dimensional parameter space that dictates the reaction kinetics, network structures, and physical properties of the final material. Therefore, it is important to have a fast way to map the polymerisation parameters to the final topological properties of the polymer network.
We have matched various idioms present in polymer chemistry to corresponding graph-theoretical analogues and indicated how these can be predicted knowing the input parameters of the system by means of analytical expressions. For instance, we gave analytical quantifications of the gelation time, the topological phase transition and the associated to it molecular weight singularity, the molecular size distribution and its asymptotes, the gyration radii of polymers, and the scaling of the monomer neighbourhood size. Some of these findings also provide an unexpected qualitative insight on chemistry of polymerisation. For instance, we have revealed the existence of two asymptotical modes that appear in the molecular size distribution, the fact that these size distributions might feature a peculiar oscillating behaviour and that the gyration radii in step-growth polymerised molecules are not dependent on time but only on the sizes of these molecules.
On a broader scope, this work paves the way to viewing polymer networks, as well as other types of network formations in condensed matter, as being complex networks in which the topology and additional layers of information (multiplexity) all contribute to the formation and function of the macroscopic material. The theoretical framework that addresses multiplexity is a topic of current network science, and one of the goals of the authors is to advance the understanding of multiplex networks in condensed matter in their future works.
Although they were produced to aid polymer chemistry in the first place, these findings are also relevant to a broader network science community as one can view polymerisation as a process that is related to percolation. In this way, the paper amounts to understanding percolation in directed networks, which turn out to feature a richer behaviour then in undirected networks. Moreover, maintaining the link to polymers and polymerisation allows one to validate complex network theories with physics. Newly developed theories of such kind have been predominantly compared to experimental data derived from single point observations, which are hard to reproduce. For example, we cannot grow the Internet from scratch again. Yet, the degree distribution may be measured by NMR and more experimental techniques that will extract network properties from matter are on their way.
We finalise the paper with a note of caution that is addressed to the whole modelling community of polymer networks. The network analysis of the node-neighbourhood scaling in the configuration model points out the existence of unphysical features that appear after the gel transition. Importantly, the unphysical scaling is not an artefact of our approach but rather an implication of the commonly trusted assumption of chemical systems being well-mixed. This assumption is standard for many modelling methods that do not account for spacial embedding of the chemical species in the three-dimensional space, as for instance is the case for the rate equations, Flory-Stockmayer theory, population balance equations, kinetic Monte Carlo, and other methods. We therefore would like to encourage the search of new network models that bring together mutual interaction of the topology and space.

Master equation of the degree distribution.
In this subsection we briefly summarise the theory from ref. 29 that allows us to recover the time evolution of the bivariate degree distribution by constructing an analytically solvable master equation. We distinguish monomer species by counting the numbers of functional groups of both types I, J and the numbers of in-and out-edges i, j. During the progress of polymerisation the functional groups are converted into chemical bonds between the monomers, and the concentration profiles M i,j,I,J (t) obey the following master equation: The master eq. (10) is a linear differential-difference equation that can be transformed to an analytically solvable system of partial differential equations by applying the generating function (GF) transform. We thus directly proceed by writing the expression for the degree distribution:  Below we make use of the expressions for μ i,j (t) to determine the gelation conversion, whereas u(i, j, t) is linked to the distribution of molecular weights, to the average molecular weight, and to the typical shortest path length. Generating functions of the degree and excess degree distributions. In this section, we omit the time dependance in the degree distribution (11), and demonstrate how one can study directed networks defined by the bivariate degree distribution as combinatorial species 55,56 . The GF of a bivariate distribution is formally given by 57 : . The excess distributions u in (i, j) and u out (i, j), are defined as the degree distributions of nodes that are reached by randomly choosing an in-or out-edge: The GFs of the latter distributions are given by: Having these expressions in hand allows us to link the moments of the degree distribution u(i, j, t) as defined by eq. (14) to the partial derivatives of the GFs by writing: These relations are used below to derive analytical expressions for various global features of the polymer network.

Molecular weight distribution.
From chemical point of view any cluster of monomers that are connected together by means of covalent bonds is considered to be a molecule. In our directed network a molecule is therefore represented by a connected component, whereas the molecular weight is simply the size of this component. The distribution of molecular weights is a popular descriptor of polymer materials. In fact, there are two ways to define such distribution: the probability w(s) that a randomly chosen monomer belongs to a component of size s is called the molecular weight distribution. Alternatively, by applying the weight s 1 we obtain the molecular size distribution, that is the probability that a randomly chosen molecule has size s. In the latter equation C provides the appropriate normalisation of probability.
Here we link the molecular weight distribution w(s) to the size distribution of connected components in the directed configuration model as derived in ref. 43 , and briefly discuss the insights that this interpretation brings to understanding the step-growth polymerisation polymerisation process. The first values of w(s) can be found by following simple considerations. For instance w(1) is the probability to choose an isolated node with no neighbours, and therefore: Furthermore, w(2) is the probability that a randomly chosen node has one edge and its only neighbour has no edges except the one that connects it with the first node: Continuing this list would lead to a combinatorial explosion of possibilities. A much faster alternative is to employ the GFs. The GF for w(s) is formally defined as: s s and is obtained from the following system of functional equations:  , and the bivariate convolution is defined as ∑ * = .
i i i j j j , 1 1 2 2 1 2 1 2 In practice, numerical values of the convolution can be conveniently obtained by Fast Fourier Transform FFT. The asymptotical analysis of w(s) yields two distinct asymptotes: The exact expression for the coefficients C 1 and C 1′ are given in ref. 43 and ν nm are as defined in eq. (14). operations, explicit analytical relations can be obtained in some special cases. Consider the case when there is only one monomer species bearing I groups of type A and J groups of type B, that is the A I B J monomer. We will now derive an explicit analytical expression for w(s). Note that in this case, the distribution of initial functionalities is trivial, P(I, J) = 1, and therefore the expression for the degree distribution given in eq. (11) simplifies to a bivariate binomial distribution:  Since u(i, j) has a binomial form, one can analytically solve the convolution powers appearing in eq. (24) to obtain:

Systems with a single monomer type.
where factors A, B, and C are defined via the Hypergeometric function: This result is identical, up to the normalisation factor s, to the size distribution derived in ref. 11 . Note we normalise with s to guarantee that w(s) is a proper probability mass function, that is ∑ = ≥ w s ( ) 1   After gel transition, the latter expression becomes more complex and reads: in in out in out o ut out in Phase transition and gel fraction. During the evolution of the network the functional groups are converted into edges and at some critical point the system accumulates so many edges that it percolates. This critical moment can be identified by a few alternative methods. For instance, one may study the asymptotical behaviour of the size distribution of connected components as given by (27). This asymptote becomes scale-free at the critical point. The other alternative is to directly detect the percolation phase transition by looking at the degree distribution itself. In this case, the changes that occur in the degree distribution at the critical point are more subtle, yet they can be detected by a specially designed criticality criterion. This criterion was introduced by Molloy and Reed for undirected networks 47 , and was later generalised to the case of directed networks in ref. 29 . In this section we briefly discuss the implications of the latter theory on our dynamic polymer network. If the only available information about a system is its degree distribution, we can detect whether the system is in the gel regime by the following criterion: The conversion of A-groups at the critical point is given by: Thus, if p A (t) > p A,critical the system contains gel. Some system, however, never produce gel. This happens because the initial configuration of the system does not have a sufficient amount of high functional monomers, or too many monomers of functionality one that terminate the growth of the network. In either case this statement can be quantified by looking at the moments of the functionality distribution P(I, J): the phase transition occurs in finite time if at least one of the following conditions is true: The gel fraction is defined as the probability that a randomly selected node belongs to the gel molecule. The GF of the component size distribution W(x) only describes the components of finite size, and the gel fraction is found as the mass deficit that departs from zero at the phase transition. The amount of this 'lost' mass, that is the probability that a randomly chosen monomer belongs to the gel, is given by f where r = W (1). This means that in order to recover g f , one needs to solve the equation for W(x) only at a single point x = 1. By substituting x = 1 into (22) one obtains: where U(x, y), U in (x, y) and U out (x, y) are given by eq. (19). It is important to note that the directionality of the network, as it is introduced in the present paper, is decisive only in the case of an asymmetric degree distribution ≠ u i j u j i ( , ) (, ). In case of a symmetric distribution u(i, j) = u(j, i), the directed network model will give the same results as the undirected model supplied with degree distribution = ∑ = + u k ui j ( ) ( , ) k i j . This implies that an AB-system (undergoing AB-reactions) with symmetric initial input (P(I, J) = P(I, J)) leads to the same network as an AA-system (undergoing AA-reactions) with = ∑ = + P K P I J ( ) ( , ) For example, the AB-system A 3 :B 3 = 1:1 forms a topologically identical polymer network (when disregarding the directionality of the edges) as the AA-system A 3 . With this in mind, eq. (38) simplifies for the case of symmetric functionality distribution P(I, J) = P(I, J) as follows: 1. Symmetry leading to with ν 01 = ν 10 and ν 02 = ν 20 . 2. Mapping the symmetric AB-system to an AA-system (e.g. A 3 :    By using induction we arrive at the following expression for the expected number of the m th -degree neighbours in terms of the degree distribution moments: The latter equality transforms to: T l T 1 1 We will now perform an asymptotic analysis of the latter equation by assuming that  l 1, in which case, the asymptotic behaviour of A l is driven by the leading eigenvalue of A. Using the eigendecomposition A l = PD l P −1 gives − = + − Note that the gel criterion given by eq. (37) can also be rewritten as a determinant: Thus, at the gel point, the matrix A′ has at least one eigenvalue equal to zero. The relation of the eigenvalues of matrix A′, λ′ 1 and λ′ 2 , to the eigenvalues of A, λ 1 and λ 2 , is as follows: λ λ ′ = − 1 1 1 and λ λ ′ = − 1 2 2 . Furthermore, above the gel transition, λ λ ′ = ′ ′ < A det 0 1 2 . This is the case, only if one eigenvalue, λ′ > 0 1 , is positive and the other one, λ′ < 0 2 , is negative, and therefore, the eigenvalues of A satisfy λ 1 > 1, λ 2 < 1 and |λ 1 | > |λ 2 |. The implications for eq. (61) are as follows: for large l, λ λ  l l 1 2 , and consequently, the total number of nodes N contained within a topological ball of radius l features the exponential growth after the gel transition:

Data Availability
No datasets were generated or analysed during the current study.