Introduction

Although it has been long hypothesised that structure of a covalently-bonded polymer network defines macroscopic physical properties of the material1, little is known about what structural patterns are being formed in such networks and what parameters decide their fate. Polymerisation is driven by reaction kinetics and physical chemistry, and at the same time, the kinetics is itself mediated by the evolution of the polymer network and its emergent geometry. Such intertwined feedback mechanisms constitute a complex paradigm for modellers. One of the classical results that dates back to the works of Flory and Stockmayer2,3,4 states that the precursor functionality has a strong impact on the network structure and therefore on its properties. However, it remains unclear to what extent is there a freedom to manipulate a network structure when the precursor functionality is fixed.

In this work, we use molecular simulations to perform a screening study of a range of monomers of identical functionality by varying the linker length and adding non-functional groups. Namely, we study difunctional (meth)acrylates, see Fig. 1, that bear two double bonds allowing them to appear in the final network as nodes with 0, 1, 2, and 4 connections. Hence the only difference between the monomers lies in the length of the linker between the double bonds and/or the presence or absence of a bulky non-functional group close to the radical site (methyl substituent). Standard theories that coarse-grain repeat units by ignoring their structural formulas postulate emergence of universal network properties for all the monomers5,6,7,8. Although such coarse-grained theories are irreplaceable when studying the rates of multiple competing reactions or tuning species concentrations in a recipe, the purpose of this work is to put this hypothesis under scrutiny and investigate if subtle changes to non-reactive parts of monomers may nevertheless reflect on the structure of the network, and therefore, on the thermo-mechanical properties of the polymer.

Fig. 1: Radical polymerisation of diacrylate monomers.
figure 1

The network grows only at a few locations called monomer radicals (labelled with a red square). When reacting, a monomer radical and another node become connected. A node may obtain maximum 4 connections as it is identified with one of the following acrylate monomers: (1) 1,10-decanediol diacrylate (DDDA), (2)1,8-octanediol diacrylate (ODDA), (3) 1,6-hexanediol diacrylate (HDDA), (4) 1,6-hexanediol dimethylacrylate (HDDMA), (5) 1,4-butyldiol diacrylate (BDDA), (6) 1,4-butyldiol dimethylacrylate (BDDMA), (7) 1,2-ethylenediol diacrylate (EDDA), (8) 1,2-ethylenediol dimethylacrylate (EDDMA).

Full size image

Geometrical and topological aspects of polymer networks can be simulated with molecular dynamics (MD), therefore enabling studies of their effects on the thermo-mechanical properties9,10,11,12,13,14. Jung et al.13 used MD simulations to study the chain-growth polymerisation of the monofunctional methyl methacrylate (MMA) and reported volume shrinkage and glass transition temperatures. Demir et al.12 and Huang et al.11 employed MD to simulate copolymerisation of the vinyl ester/styrene polymer network. These authors reported structural polymer properties, such as the formation of cycles, the glass transition and Young’s modulus. Klähn et al.10 studied the effect of different monofunctional (meth)acrylate monomers on the glass transition temperature using molecular simulations. Rudyak et al.9 studied the effect of the initiator concentration on a chain-growth polymer network and reported gel conversions, gel fractions and cycle length distributions. In our previous study, we established a simulation protocol for polymerisation of 1,6-hexanediol diacrylate (HDDA), and observed the formation of ‘short cycles’ during the polymerisation process: small cycles were created in primary and secondary cyclization reactions and were promoted by enhanced monomer flexibility15. On the one hand, the ability to form cycles is dictated by the local conformational dynamics of monomers, namely by the torsional strain and steric hindrance. On the other hand, these reactions influence the global properties of the network, such as the density and excluded volume by introducing topological defects that are elastically inactive. This has motivated us to enquire whether chemically different monomers with identical number of functional groups may lead to different gelation times and elastomeric properties of the network. Since the non-functional part of monomer precursors can be chemically manipulated, the mechanism linking the monomer structure with the final polymer network is important for enabling, for example, the rational design of materials.

The problem of optimising molecular weight, elastic modulus, or glass transition temperature arises in different contexts, for example, when designing self-healing materials16, but also in lithography, coatings for biomaterials, textiles and dental composites17,18,19,20,21,22,23. Polymer resin recipes are designed to include multiple multifunctional monomer precursors to achieve a desired degree distribution8. Acrylate polymers assembling into covalently-bonded crosslinked networks with several vinyl groups per monomer24 and can be tuned by adjusting different substituents in the α-carbon, linkers between the vinyl groups, and by adding functional groups25. In addition, one can manipulate chemical kinetics by adding a chain transfer or chain shuttling agents to reduce molecular weight24, and these techniques may also be expected to affect the overall network structure.

In the Results section, we report physical properties related to monomer diffusion, glass transition and elasticity. Then, we turn to the topological analysis of the network and its emergent geometry and match the structural changes in the network with the development of the physical properties. We show universality across the studied monomer types with respect to such properties as the degree distribution, cyclomatic number, gel point conversion, and Young’s modulus. We also report that such universality is broken in the dense networks, as can be seen by diverging glass transition temperatures associated with different monomer types. We show that this divergence may be mediated by small chordless cycles, also called topological holes, and report two qualitative transitions in the network structure. The first transition—gelation, is marked by the appearance of an extensive cluster and is responsible for emergence of elastic behaviour. The second transition is marked by emergence of an extensive higher-order topological entity called a cell complex, and is associated with a steep increase of glass transition temperature that may cause an onset of brittleness.

Results

Polymerisation process

During radical polymerisation, initially disconnected monomer units join together to form a growing network through three main reactions: initiation, propagation and termination. We schematically visualise these processes by representing monomers as labelled nodes of various functionality that receive connections according to the rules shown in Fig. 1. The real radical polymerisation is a process that is also affected by the microscopic properties, e.g. monomer mobility and geometry, as well as by the macroscopic topological and thermo-mechanical changes in the whole system, e.g. viscosity, gelation, glass transition, volume and overall density. We study this process with dynamical systems that govern the motion of all atoms individually. The methods section explains the set-up of our curing simulations and the validation of the force field model.

Diffusion

In growing polymer networks, reactions rapidly become diffusion limited. Bond conversion 0 ≤ χ(t) ≤ 1 is the ratio between the formed covalent bonds at time t and the maximal number of covalent bonds that the system may contain. Since each monomer can have maximum four neighbours, and each radical reduces the total number of bonds by 2, we have:

$$\chi (t)=\frac{b(t)}{2{N}_{\text{sys}}-{R}_{0}},$$
(1)

where b(t) is the current number of bonds, R0 = 100 the initial number of radicals, and Nsys = 2000 is the number of monomers in the system. Figure 2a quantifies the kinetic slowdown χ(t) for each monomer type, where we observe that: (1) for a given linker length, the presence of a methyl substituent significantly slows down bond formation, (2) for a given substituent, the polymerisation process slows down with increasing linker length, and (3) the mobility of all systems slows down with curing. In Supplementary Fig. 4, we compare the diffusion coefficients of free monomers as a function of conversion for systems with and without a methyl substituent. These results show that the decrease in the speed of polymerisation is affected by the overall reduced translational mobility of bulkier dimethacrylates. One can also conclude from the radial distribution functions (RDFs) in the melt state, shown in Supplementary Fig. 5, that the addition of a methyl substituent reduces the frequency of reactions, as reactions involving the α-carbon become less likely due to the steric effects. The analysis of the RDFs suggests that shielding a radical site by longer linkers is negligible, see Supplementary Fig. 6, which was also pointed out by Kurdikar and Peppas26. Therefore the diffusion of free monomers might play a dominant contribution to the overall kinetic slowdown, see Supplementary Fig. 7, which constitutes perhaps the strongest difference between the studied monomers.

Fig. 2: Inferring free monomer diffusion and thermo-mechanical properties from molecular simulations.
figure 2

a Bond conversion vs. simulation time. Confidence intervals indicate one standard deviation. Inset shows time in logarithmic scale. b The procedure for calculating Tg as the break point in a piecewise linear fit of density vs. temperature, see Methods for details. c Tg vs. bond conversion features a universal collapse to a linear dependence. Supplementary Fig. 14 reports confidence intervals for these data. d The stress-strain curve with a thermal noise is processed to identify the region of initial linear growth and its slope, see Methods for details. The deviations from the fit, shown in the inset, are confirmed to be normally distributed and not autocorrelated. e Young’s modulus features a linear dependence on bond conversion resulting in fuzzy linear fits due to thermal fluctuations. The cross-sections of the probability densities of the fits at χ = 0.8 are shown in the inset. Supplementary Fig. 15 reports confidence intervals. f Prediction of Young’s modulus from solely spectral graph properties of the network, is shown for the HDDA monomer. See Supplementary Fig. 17 for resistance-elastic moduli comparison for other monomers. The elastic modulus and Tg were averaged over four simulation runs.

Full size image

The polymerisation rate of methacrylates was observed to be lower than that in acrylates for mono- and high-functional monomers, both in experiments27,28 and computationally29. In previous studies on monofunctional acrylates, the decrease of the propagation rates with increasing linker length is often attributed to the increase of the activation energy. As we only observe effects of diffusion and steric hindrance, it is more accurate to compare the trend in Fig. 2a with reported trends of the pre-exponential factors of the rate coefficients. Mavroudakis30 reports the pre-exponential factors for monomers with varying linker lengths, which is in good agreement with our findings. Kurdikar and Peppas26 investigated the effect of the linker length on the reaction kinetics of the polymerisation of different poly (ethylene glycol) diacrylates and found a decrease in the rate constant for longer linkers, attributing this effect to a decreased monomer diffusivity.

After confirming that chemical reactions take place on different time scales for different monomers, in what follows we report all results as a function of bond conversion χ instead of time. In this way we focus on the outcome of the polymerisation and not merely on its absolute speed.

Thermomechanics

Tg is determined by analysing the density as a function of temperature31. As illustrated in Fig. 2b and explained in Methods, we search for such a temperature, where this function fails to be smooth. As shown in Fig. 2c, Tg increases with bond conversion for all monomers. The increase is first linear and well-explained by a universal relationship across all monomer species:

$${T}_{\text{g}}(\chi )=A\chi +{T}_{\text{g,0}},\ A=155\,\text{K}\,,$$
(2)

when χ < 0.7. The linear trend is abruptly broken at higher bond conversions. In the Topology section, we argue that a novel higher-order percolation transition provides an explanation for such an abrupt change of Tg. The data presented in Fig. 2c suggest that shorter linkers might promote higher Tg for χ > 0.7.

Both observations, the linear dependence and the deviation that follows it, are in line with findings in literature. Bowman et al.32 reported the Tg for the homo- and copolymerisation of DEGDA and DEGDMA with monofunctional monomers (n-octyl methacrylate and n-heptyl acrylate, respectively). These authors observed an increase in Tg with increasing crosslinker density. Bowman et al.33 also studied the effect of the distance between functional groups for poly (ethylene glycol) dimethacrylates from 1 to 4 ethylene glycol units and reported glass transition between 250 and 280 C. Kurdikar and Peppas34, reported the Tg for the different poly (ethylene glycol) systems and reported that a longer linker reduces the Tg. Klän et al.10 also observed lower Tg for longer linkers in polymerised (meth)acrylates using molecular simulations. Jerolimov et al.35 studied the effect of divinyl crosslinking agents with different chain length (EGDMA and TEGDMA) on the Tg of PMMA and also reported a similar observation. Recently, Klän et al.10 observed that Tg increases in the presence of the methyl group using molecular simulations.

Elasticity

Young’s modulus is determined by estimating the slope of the stress-strain curve in the linear region during of the static tensile test. We performed the test under conditions of thermal fluctuations, see Fig. 2d, and statistical tests were applied to infer the linear region and its slope (see Methods).

Figure 2e shows that for all systems Young’s moduli increase linearly as a function of bond conversion:

$$E(\chi )=B(\chi -{\chi }_{0}),$$
(3)

featuring a universal slope B = 3.6GPa, with χ0 being the largest χ, for which E = 0. Parameter χ0 turns out to be weakly dependant on linker length: Monomers with longer linker length feature smaller values of the moduli with an exception of HDDA-BDDA pair that deviates from the universal behaviour for χ < 0.6 showing a reverse order. The overall effects of the crosslinker length and bond conversion on elasticity are in agreement with experimental findings36,37.

The evolution of the elastic modulus is universally explained for all monomers also by a spectral property of the networks called graph resistance, see Methods for calculations. Figure 2f compares the modulus obtained from the tensile test and the scaled graph resistance, for example of HDDA. As shown in Supplementary Fig. 17, presenting this comparison for other monomers, the graph resistance has to be scaled by a constant that absorbs the contribution from a single monomer coil and ranges between 5.0 and 8.8 GPa.

Topology

As a result of polymerisation, monomers join together with covalent bonds forming dense polymer networks, for example, as the series of snapshots shown in Fig. 3a for the DDDA system. So far, we have discussed emergent physical properties of the polymer related to monomer diffusion, density and elasticity. For all of these quantities, with an exception of the monomer diffusion coefficient, the number of bonds was a more important determining factor than the chemical type of the monomer. This observation suggests that the physical properties under investigation are mainly defined by the network structure, and only indirectly, by the monomer type. In other words: the network is a necessary intermediary link between the monomer structure and the macroscopic physical properties.

Fig. 3: Two-dimensional layouts of networks and cell-complexes developing in polymerising DDDA.
figure 3

a Network structures with highlighted components. The first percolation transition is thought to take place at χ = 0.2 and shortly afterwards the whole system predominantly consists of one connected component (isolated nodes not shown). The largest component continues to evolve as new bonds appear internally. b Evolution of the largest connected component in the cell complex. Nodes represent holes in the former network: node size represents the size of the hole, with smallest nodes corresponding to size 3 and the largest to size 6. The second percolation transition is thought to take place at χ = 0.7, wherein the whole cell complex becomes spanned by a component of an extensive size, which severely restricts mobility of the polymer network.

Full size image

We quantify the structure of polymer networks in terms of several notions used in topology and (spectral) graph theory. Namely, we look at degree distribution, graph resistance, connected components, cycles, holes, and cell complexes. We analyse the dependence of these quantities on the monomer type, and argue about their effect on the emergent physical properties.

Our first finding is that self loops form frequently, as caused by a radical reacting with the pending vinyl group of the same monomer. Each self-loop decreases the maximum degree of one monomer from 4 to 2, and therefore self loops are of general interest in reaction kinetics and material science, where they are often referred to as topological defects38,39,40,41,42,43. In a similar fashion, a double edge occurs when a pair of monomers reacts twice, thus becoming connected with two bonds, which limits their external number of neighbours to 4 per pair, instead of 6. As Fig. 4a suggests, a shorter linker between the functional groups promotes formation of self loops and a very similar trend is observed for double edges, see Fig. 4b. This observation holds with one exception of HDDA. The ab initio and classical simulations presented in Supplementary Figs. 2 and 13 suggest that this exception might be related to a coiled metastable state of HDDA—the end-to-end distance being close to the reaction cutoff radius. Moreover, the ability to form triangles—cycles of size three—is also enhanced by short linkers, see Supplementary Fig. 10, again with the BDDA-HDDA pair being swapped in the order of the trend.

Fig. 4: Graph-theoretical and topological properties of the polymer network.
figure 4

a Number of nodes with one self-loop. b Number of nodes with either a double edge or a self loop; the inset shows an averaged distribution of hole sizes at χ = 0.8. c Universal collapse of the overall cyclomatic complexity. d The fraction of holes of indicated size as a function of bond conversion. Supplementary Fig. 12 reports hole statistics for each monomer separately. e Degree distributions as a function of conversion feature universality. f Universal collapse of the sizes of the largest component and cell complex in different systems. The largest connected components traverses percolation transitions between χ = 0.15 and χ = 0.25 depending on the monomer type. The sizes of the largest connected cell complex feature a second percolation transition around χ = 0.7. Supplementary Fig. 11 reports the confidence intervals.

Full size image

As shown in Supplementary Fig. 9 the effect of the methyl substituent is less clear: For EDDMA, the methyl substitution induces more loops, but for the other two linker lengths (BDDMA and HDDMA) the results are inconclusive. The end-to-end distance probability distribution of the pure monomers, shown in Supplementary Fig. 2, suggests that the larger number of self loops in EDDMA could also be caused by a metastable coiled configuration of this molecule.

Such a strong influence of monomer types on forming cyclic structures is characteristic only to the smallest scale. To demonstrate this, we compute the maximum number of edges one can remove before the network becomes composed of tree components—the cyclomatic complexity r, which can be calculated from the relationship r = M − N + C with M being the total number of edges, N – number of nodes, and C – number of connected components in the network. Clearly, self loops, double edges, and triangles increase cyclomatic complexity. As shown in Fig. 4c, all monomers lead to a comparable cyclomatic complexity, where about 10% of this quantity is explained by the self loops and double edges alone. In order to investigate higher order cyclic structures, we must differentiate between a cycle, which is a closed path in a graph, and a hole—a cycle for which any subpath is also the shortest path connecting its own endpoints.

Bonds that are not self loops, connect different nodes together and thus extend the network in its size. By computing the degree distribution, that is fractions of nodes having different numbers of neighbours, we again observe a universal behaviour across monomer types, shown in Fig. 4e. The dynamics seen in Fig. 4e is different from the binomial degree distribution, which is characteristic to the standard bond percolation5,7,8,44, and resembles the chain-growth polymerisation45. The degrees of monomers with methyl substituents are analysed in Supplementary Fig. 8, showing that the presence/absence of the methyl group does not have a significant effect on the degree distribution. This is in strong contrast to the observed differences in polymerisation properties for different monomer types when studied as a function of time, for example as in Fig. 2a, where such differences are mainly attributed to diffusion limited nature of the reaction. Figure 4e gives reason to suggest that the diffusion effects do not strongly influence the density of the network when viewed as a function of bond conversion.

Self loops and small cycles are in competition with forming bonds that expand the network size. Although formation of small cycles does not significantly affect the degree distribution, they delay the onset of gelation, that is the moment in time when an extensive connected component emerges. Note that the mean-field theories predict the gel point exclusively on the basis of the degree distribution5,7,8,46,47, which unlike self loops and small cycles, was shown to be insensitive to the monomer type. The onset of gelation is identified by monitoring the fraction of nodes in the largest connected component. As shown in Fig. 4f, the gelation occurs later in bond conversion for shorter linkers. This trend is attributed to a large number of self loops that act as ‘defects’—they increase the bond conversion without contributing to the expansion of the network. A similar conclusion about the effect of cyclization on the gel point was reported by Elliott et al.48,49, where the gel point of TeGDMA was observed to be higher than the gel point of BisGDMA (bisphenol A-glycidyl dimethacrylate), with TeGDMA exhibiting threefolded more cyclization than BisGDMA due to the differences in flexibility of the monomers.

Since the gelation starts rather early, around χ = 0.2, most of the network evolution takes place with the giant component being present. Moreover, around χ = 0.4, the whole network predominantly consists of one connected component, plus isolated nodes. This behaviour is illustrated with network snapshots in Fig. 3a. Surprisingly, the connected network undergoes another structural transition, long after the onset of gelation. We demonstrate this by analysing the cell complex50,51,52, which builds upon the idea of a topological hole. Our construction of the cell complex can be thought of as a higher order network in which the nodes are identified with the holes in the original polymer network, and a pair of such cell-nodes is connected if the corresponding holes share vertices53, see Fig. 3b. This definition allows us to identify connected components in the cell complex and calculate their size. The sizes of the largest component in the cell complex, see Fig. 4f, suggest that such a construction also undergoes a percolation transition, around χ = 0.7, when a percolating cell component emerges. By χ = 0.8 − 0.9, the giant cell complex expands further and incorporates the majority of the nodes. The structure and development of the cell complex resembles a tree-like construction with a small number of higher order cycles, which can be compared to the early stages of the network formation. There is a large difference between a monomer-node in the original network and the hole-node in the cell complex—because small holes limit the uniaxial rotation of the constituting monomers—the latter has less degrees of freedom. A rapid change in the degrees of freedom that occurs in the physical network at the second transition, see Fig. 5, may provide a possible reason behind a sudden nonlinearity in Tg(χ) observed in Fig. 2c for all monomer systems. The location of the second transition is affected by sizes of small topological holes, which, in turn, are influenced by monomer conformation dynamics. This mechanism allows monomers of different types to affect Tg, albeit only at relatively late bond conversions after the onset of the second transition.

Fig. 5: Mean square displacement (MSD) of atoms in the EDDA system during curing.
figure 5

The figure demonstrates a separation of the network into two phases with respect to thermal fluctuations. On average, the MSD of all atoms decreases with curing, however, the atoms that eventually join the largest connected component in the cell complex (identified at χ = 0.6) reach the terminal MSD earlier than the average atom in the system.

Full size image

Discussion

In this work we have re-evaluated the causal link between the non-functional form of the monomer precursor and the physical properties of the resulting polymer material. To this end, we studied diacrylate monomers with different linker length and a present/absent bulky methyl group to show that in the causal chain between monomer structure and the emergent physical properties, the network topology plays a role of a unavoidable intermediary.

Monomer preference to occupy certain coiled states is central to such a mechanism as it affects formation of topological holes. We do not entirely confirm the random-coil trend54, longer linker—less loops. For instance, despite being longer than BDDA, the HDDA monomer has a stronger preference to form loops, as was shown by classical and also confirmed by the ab initio metastability simulations.

Although the degree distribution is only weakly affected by the monomer type, the way the edges are arranged in the network is strongly affected as can be seen from the size distribution of topological holes. These differences affect the number of edges it takes to reach the gelation—the moment when an extensive connected component appears and transforms formerly viscous liquid material into an elastic solid. We observed such a transition in two ways: topologically—by monitoring the largest connected component in the network, and physically—by detecting divergence of the tensile modulus from zero. Even though longer linkers were shown to slightly accelerate gelation, the tensile modulus was found to be smaller for monomers with longer linkers.

We showed that although Tg is not affected by gelation, there exists a novel higher-order percolation transition that drastically ramps up the Tg at late bond conversions. This percolation transition marks the appearance of an extensive cell complex that spans the network and therefore constrains its mobility. In a similar way as gelation marks the onset of elastic behaviour in the polymer, the higher order percolation transition makes the material more brittle, as was shown by the Tg tests. Since Tg suddenly becomes sensitive to monomer type at late stages, we postulate that one should be capable of influencing glass transition temperature by manipulating non-functional parts a monomer, e.g. by changing linker length or adding non-reactive groups. A parallel can be made here with another type of higher order structures constraining mobility of physical networks—a k-core—that is related to the onset of jamming in granular gels55.

Even though we observe strong differences in the polymerisation kinetics as a function of time for different monomers, the topological and physical properties are surprisingly similar when reported as a function of conversion. Universal behaviour was seen in almost identical degree distributions, numbers of large cycles, and sizes of components far from the percolation transition. Both Tg and Young’s modulus feature linear dependence on the bond conversion with a universal slope at early to intermediate stages of polymerisation. Therefore, depending on what quantities of interest one aims to predict, coarse-grained and macroscopic models may be well justified, for example, when used in the partial differential equations absorbing MD generated data56,57.

Methods

MD model and system set-up

The monomers were modelled at the united-atom (UA) level with the TraPPE-UA force field by Maerzke et al.58 developed for acrylates. The covalent bonds were initially represented by a harmonic potential as opposed to rigid bonds, which is necessary for the crosslinking simulations. See Tables 1–4 in the Supportive Information in reference15 for the functional form of the force field.

The initial systems were set-up by randomly packing 2000 monomers in a 3D simulation cell (initial size 200 × 200 Å) with periodic boundary conditions. The systems were then equilibrated in the NVT ensemble at 600 K for 10 ps using a time step of 0.1fs and further equilibrated in the NPT ensemble at 600 K for 500 ps with a time step of 1fs. This resulted in three-dimensional periodic simulation boxes of around 100 × 100 Å. The simulations were carried out using the Nosé–Hoover thermostat59,60, and, in the constant pressure cases, using the Martyna–Tuckerman barostat61. From the 2000 monomers (which amounts to, e.g. 32000 HDDA atoms) units, 5% were considered active at the beginning.

To validate force field and study the properties of the liquid monomers at ambient conditions, such as density, self-diffusivity and the viscosity of the system, the equilibrated systems at 600 K were cooled down using the NPT ensemble from 600 K to 300 K and further equilibrated for 500 ps at 300 K. The equilibrated systems at 300 K were then used as input for longer simulations at 300 K in the NVT ensemble to determine the density, self-diffusivity and viscosity properties of the melts, which are in a good agreement with experimental values, see Supplementary Table 1. For the viscosity, a correlation time of 1 ns was used and the simulations were run for 500 ns. Further validation of the force field was done by comparing histograms and RDFs of HDDA to ab-initio based MD as implemented in the CP2K software62. The ab initio simulations were performed using the BLYP functional supplemented with D3 dispersion corrections25,63. As shown in Supplementary Figs. 1 and 2, the histogram of the dihedral angles of HDDA in liquid state and the end-to-end distance in HDDA are well-captured by the TraPPE-UA force field with a few notable differences. In addition, the RDFs of the 1-1 and 2-2 reactive sites are also well-captured by the employed force filed in comparison to the ab initio simulations as shown in Supplementary Fig. 3.

Curing simulations

We used a simulation protocol that accelerates the curing process by lowering activation energies and increasing temperature to 600 K, which was earlier shown not to affect the structure of the network when considered as a function of bond conversion15. The covalent bonds appearing during the curing of the network were created using a reactions cutoff radius of 4 Å. The proximity of reactive sites was checked every 10 fs and a bond was created with probability 0.5 for every time reactive sites within the cutoff radius. After placing a new crosslink, the whole system was equilibrated in the NPT ensemble. The activation energy, which affects the chemical reaction rate, was assumed to be equal for all monomers. Hence, to link the simulation time to the real time, one needs additionally to compensate the contribution of the chemical activation energy. The curing was simulated for up to 100 ns, which is enough for all systems to reach bond conversion of 0.8 in our protocol. This corresponds to 4 weeks of computational time for one trajectory. The curing simulations and the tests for physical properties were repeated four times to average out the effect of the initial condition bias. To validate the protocol, we computed the volume shrinkage factor during curing, see Supplementary Fig. 18, which matches experimentally observed values.

The curing of (meth)acrylates involves the addition of radicals to the unsaturated double bonds in the vinyl groups. Since the initiation reaction is much more slower than radical propagation, we assume that a 5% of radical sites is active from the start. In experimental conditions radicals are either continuously initiated or appear in bursts as controlled by the pulses of light irradiation. This results in a higher radical concentration than typically used in experiments, however lowering the radical concentrations would make simulations impossibly long. For the effect on the initial amount of radicals on the network formation, see Section 4 in ref. 15 and Figs. S6 and S7 therein. Assuming the initial reactive radical sites are specified, the curing was implemented in three consecutive steps: (1) the addition of a bond to the unsaturated carbon that undergoes a reaction, (2) the change of the carbon character from unsaturated to saturated (to prevent it from participating in future reactions), and (3) the propagation/regeneration of the reactive site by changing the character of the neighbouring carbon from unsaturated to the reactive one.

Calculations for T g

The glass transition temperature was obtained by performing NPT simulations for a series of different temperatures. The systems were cooled down from 600 K to 150 K at the rate of 5 ⋅ 109 K s−1. We determined Tg by analysing the density as a function of temperature31 and performing a piecewise linear fit. In order to identify Tg and quantify the uncertainty that arises due to the finite system size and the thermal fluctuations we introduce the following protocol for statistical inference: First, we solved the minimisation problem

$${j}^{* }=\arg \mathop{\min }\limits_{j}\mathop{\sum }\limits_{i=1}^{j}\parallel {r}_{i,1}{\parallel }_{2}+\mathop{\sum }\limits_{i=j+1}^{n}\parallel {r}_{i,2}{\parallel }_{2}$$
(4)

where

$${r}_{i,1}= \, {\rho }_{i}-({\beta }_{j,1}{T}_{i}+{\varepsilon }_{j,1}),\\ {r}_{i,2}= \, {\rho }_{i}-({\beta }_{j,2}{T}_{i}+{\varepsilon }_{j,2}),$$
(5)

are the residuals of least square linear fits that correspond to the first (ρi, Ti, i = 1, …, j) and second (ρi, Ti, i = j + 1, …, n) fragments of data points separated by index j; ρi, i = 1, …, n are the measured density at temperatures Ti, and βj,k, εj,k, k = 1, 2 are the coefficients of the fits. We then define the glass transition temperature as the point where the linear fits intersect:

$${T}_{\text{g}}=\frac{{\varepsilon }_{j,2}-{\varepsilon }_{j,1}}{{\beta }_{j,1}-{\beta }_{j,2}}.$$
(6)

Even though ρi are normally distributed around their fits, Tg is defined as a quotient, and therefore, is not normally nor symmetrically distributed. In fact the distribution for Tg may even contain a heavy tail when the intersection angle is close to π. We thus estimate the confidence intervals for Tg by using a sampling method. Let us denote the averages of the first and second intervals as \({\langle {x}_{i}\rangle }_{1}:= \frac{1}{{j}^{* }}\mathop{\sum }\nolimits_{i = 1}^{{j}^{* }}{x}_{i}\) and \({\langle {x}_{i}\rangle }_{2}:= \frac{1}{n-{j}^{* }+1}\mathop{\sum }\nolimits_{i = {j}^{* }+1}^{n}{x}_{i}\). We define

$${\sigma }_{1}^{2}:=\, {\langle {T}_{i}\rangle }_{2}^{2}\frac{{\langle {r}_{i,1}^{2}\rangle }_{1}-{\langle {r}_{i,1}\rangle }_{1}^{2}}{{\langle {T}_{i}^{2}\rangle }_{1}-{j}^{* }{\langle {T}_{i}\rangle }_{1}^{2}}\\ + \,{\langle {T}_{i}\rangle }_{1}^{2}\frac{{\langle {r}_{i,1}^{2}\rangle }_{2}-{\langle {r}_{i,1}\rangle }_{2}^{2}}{{\langle {T}_{i}^{2}\rangle }_{2}-(n-{j}^{* }+1){\langle {T}_{i}\rangle }_{2}^{2}}$$
(7)

and

$${\sigma }_{2}^{2}:= \frac{{\langle {r}_{i,1}^{2}\rangle }_{1}-{\langle {r}_{i,1}\rangle }_{1}^{2}}{{\langle {T}_{i}^{2}\rangle }_{1}-{j}^{* }{\langle {T}_{i}\rangle }_{1}^{2}}+\frac{{\langle {r}_{i,1}^{2}\rangle }_{2}-{\langle {r}_{i,1}\rangle }_{2}^{2}}{{\langle {T}_{i}^{2}\rangle }_{2}-(n-{j}^{* }+1){\langle {T}_{i}\rangle }_{2}^{2}}.$$
(8)

Let \({X}_{i} \sim {\mathcal{N}}({\varepsilon }_{j,2}-{\varepsilon }_{j,1},{\sigma }_{1}^{2})\) and \({X}_{i} \sim {\mathcal{N}}({\beta }_{j,1}-{\beta }_{j,2},{\sigma }_{2}^{2})\) be N independently generated samples generated from the normal distribution \({\mathcal{N}}(\mu ,{\sigma }^{2})\) with mean μ and variance σ2. Then \({Z}_{i}=\frac{{X}_{i}}{{Y}_{i}}\) are random samples from the distribution for the Tg estimate. The confidence intervals for Tg are the percentiles of the empirical density function for Zi − 〈Zi〉.

Calculation of Young’s modulus

To determine Young’s modulus, we performed uniaxial tension loading simulations by increasing the length of the simulation cell along the loading direction at every MD step while maintaining atmospheric pressure in the transverse directions using a barostat. Young’s modulus was determined by estimating the slope of the engineering stress-strain curve in the linear region, when a sample is pulled at a constant rate in the static tensile test. The stress-strain curves were obtained by applying a uniaxial deformation along the x-axis at a constant rate and measuring the strain-stress response curves. The limit of the longitudinal strain was set to be 50%, and the strain rate was varied between 108 and 1010s−1 showing no statistically significant dependence of the inferred Young’s modulus on the rate, see Supplementary Fig. 16. The strain was determined from \(L(t)={L}_{0}\times \exp (\,\text{rate}\,\times {\rm{d}}t)\), where L(t) and L0 are respectively the instantaneous and initial length of the simulation box in the direction of elongation. The stress was calculated from the virial expansion. The transverse directions were maintained at atmospheric pressure using the Martyna–Tuckerman barostat61. Since the test was performed under conditions of thermal fluctuations, as illustrated in Fig. 2d, we used Ljung-Box Q-Test to determine if fluctuations around the linear fit for the first j = 2, 3, . . . data-points are pair-correlated. Then, if the identified fragment of the data series also passed the Kolmogorov–Smirnov test for normality, the slope and its confidence interval were deduced applying the linear regression. The linear region was identified as the largest strain for which these tests were passed. The length of the linear region typically fluctuated between 0.06 and 0.2, with higher bond conversion being associated with shorter linear regions. This procedure was repeated for different conversions χ and different monomer types, as reported in Supplementary Fig. 15. All data points together with their confidence intervals across the conversion range were fitted again with a straight line, as shown in Fig. 2e. The ‘fuzziness’ of the latter fits is due to confidence intervals of data points in the latter fits. Such secondary fits are therefore drawn from a probability distribution. Note that in the viscous regime (before the onset of gelation), our method inferred no values for the modulus because the data violated Ljung-Box Q-Test.

Average resistance model for elasticity

The approximation of the elastic modulus was obtained using the graph resistance model:

$${E}^{* }=C{\left(\frac{1}{| {{{\Omega }}}_{1}| | {{{\Omega }}}_{2}| }\mathop{\sum}\limits_{i\in {{{\Omega }}}_{1},j\in {{{\Omega }}}_{2}}({L}_{i,i}^{+}+{L}_{j,j}^{+}-{L}_{i,j}^{+}-{L}_{j,i}^{+})\right)}^{-1}\,\,,$$

where C is a constant that represents monomer tendency to coil and (was estimated from molecular simulations), L = D − A is the graph laplacian with D being the diagonal matrix with node degrees on the diagonal and A the adjacency matrix of the monomer network and L+ indicating its Moore–Penrose inverse; Ω1 and Ω2 are the index sets for monomers that are located at the two opposite faces of the simulation box. This estimate is based on the electrical engineering concept of effective resistance64. Note that the method is not applicable when there is no percolating cluster between Ω1 and Ω2.