## 2 Materials and Methods
### Dissipative Particle Dynamics Approach
To minimize unphysical topological crossings of bonded polymer chains (a known limitation of DPD), we adapted the modified Segmental Repulsive Potential (mSRP) [47; 48; 49] formulation. Notably, the mSRP DPD formulation captures the effects of entanglements in polymer melts [47]. Within the mSRP DPD framework, pseudo beads are introduced at the center of each bond for all the bonds initially present in the system; the pseudo beads interact with other pseudo beads with a soft repulsion potential defined as [47]:  
\[\mathbf{F}_{ij}^{SRP}= \begin{cases}b\big{(}1-d_{ij}/d_{c}\big{)}\mathbf{e}_{ij}^{S}& \big{(}d_{ij}<d_{c}\big{)}\\ 0&\big{(}d_{ij}\geq r_{c}\big{)}\end{cases}  \tag{2}\]  
where \(b\) defines the strength of the interaction between the pseudo beads, \(d_{ij}=\big{|}\mathbf{d}_{ij}\big{|}\) is the distance between these beads located at the centres of the bonds, \(\mathbf{e}_{ij}^{S}=\mathbf{d}_{ij}/d_{ij}\), and \(d_{c}\) is the cut-off distance for the mSRP interactions. Herein, the mSRP framework is modified [48] to create additional pseudo beads for the bonds created during crosslinking. For the mSRP potential, we set [47]\(b=80\,k_{B}T/r_{c}^{2}\) and \(d_{c}=0.8\,r_{c}\); these parameters are shown to effectively minimize the bond crossings [47; 49]. The strength of the randomforce is chosen as \(\sigma=3\), the number density in simulations is \(\rho=3/r_{c}{}^{3}\), and the time step is \(\Delta t=0.02\ \tau\), where \(\tau\) is the reduced unit of time [27]. The \(r_{c}\), temperature, and the mass of a DPD bead are set at 1.0 in reduced DPD units [27; 45]. With the above choice of the degree of coarse-graining, the dimensionless units of length and time correspond to [45]\(r_{c}\ \approx 0.65\) nm and \(\tau\approx\ 88\) ps, respectively. All the results below are provided in the reduced DPD units. PHMS is represented by a linear chain that consists of 30 DPD beads (\(N_{\rm PHMS}=30\)) and we vary the number of beads comprising the \(v\)-PDMS chains, \(N_{v-\rm PDMS}\). The size of the simulation box is \(L\times L\times L\), where \(L=40\) in reduced DPD units as provided above, periodic boundary conditions are applied in all directions. With the chosen density of 3, the total number of beads in the system is \(1.92\times 10^{5}\). Prior to the production runs, the blends are equilibrated for \(6\times 10^{6}\) simulation timesteps with the same value of the repulsion coefficient (\(a_{ij}=78\)) chosen between all the beads within each system; no reactions occur during the equilibration. At the initial timestep of the production run, the repulsion parameters are modified between the beads using the values of \(a_{ij}\) provided below. The equations of motion are integrated using the DPD module as implemented in the LAMMPS simulation package [50; 51] with mSRP code [47] and the most recent modification of the mSRP code allowing one to use mSRP approach simultaneously with breaking or creating polymer bonds [48; 49]. All simulation snapshots were visualized using the Visual Molecular Dynamics (VMD) software [52], and Wolfram Mathematica software [53] was used to generate all the plots.