Nanostars planarity modulates the rheology of DNA hydrogels

In analogy with classic rigidity problems of networks and frames, the elastic properties of hydrogels made of DNA nanostars (DNAns) are expected to strongly depend on the precise geometry of their building blocks. However, it is currently not possible to determine DNAns shape experimentally. Computational coarse-grained models that can retain the correct geometry of DNA nanostars and account for the bulk properties observed in recent experiments could provide missing insights. In this study, we perform metadynamics simulations to obtain the preferred configuration of three-armed DNA nanostars simulated with the oxDNA model. Based on these results we introduce a coarse-grained computational model of nanostars that can self assemble into complex three dimensional percolating networks. We compare two systems with different designs, in which either planar or non-planar nanostars are used. Structural and network analysis reveal completely different features for the two cases, leading to two contrasting rheological properties. The mobility of molecules is larger in the non-planar case, which is consistent with a lower viscosity measured from Green–Kubo simulations in equilibrium. To the best of our knowledge, this is the first work connecting the geometry of DNAns with the bulk rheological properties of DNA hydrogels and may inform the design of future DNA based materials.


I Details of Metadynamics simulations
Reminder of the theory As described in the main text, Metadynamics (MetaD) is a computer simulation method that can be used to estimate the free energy landscape (FEL) of a system. Importantly, it is assumed that this FEL can be written in terms of few Collective variables (CVs), ζ i (r), that depend only on the positions, (r), of particles and that can distinguish between any two (or more) states of the system. Different states are represented by minima in the FEL.
In MetaD, the sampling of the phase-space of a system is facilitated by the introduction of a history-dependent bias potential that forces the system to migrate from one minimum of the FEL to the next one. When just one CV is considered the potential, constructed as a sum of Gaussians at time t, has the form: where B, ∆ and τg are constants representing the Gaussian height, width and frequency at which the Gaussians are added, respectively. ν(t ) = ζ(r(t )) represents the the value of the CV at time t , and t < t must be satisfied. If B is large, the free energy surface will be explored at a fast pace, but the reconstructed profile will be affected by large errors. Instead, if B or τg are small the reconstruction will be accurate, but it will take a longer time.
The basic assumption of metadynamics is that after a sufficiently long time, lim t→∞ U G (ζ, t) provides an estimate of the free energy landscape, F (ζ), modulo a constant. However, since the potential of the system changes every time a Gaussian is added, U G (ζ, t) oscillates around F (ζ) at long times. Another method that ensures that the FEL converges more smoothly is the Well-Tempered Metadynamics (WTMetaD). In this variant of MetaD, the height of the Gaussians added is no longer a constant but it scales at each timestep, so the bias potential is obtained by replacing B in Eq. S1 by: where ∆T has the dimension of temperature. This parameter is related to the bias factor (γ = (T + ∆T )/T ), which indicates how fast the Gaussian height decreases in a simulation. The higher the value of γ, the slower the height decrease. Therefore, γ → ∞ represents normal MetaD simulations and γ = 1 represent ordinary molecular dynamics (MD) simulations.

The collective variable
Here we define a collective variable (dp, described in the main text) that can distinguish between planar and non-planar geometries of the DNAns. In practice, dp is obtained from the direction of each of the dsDNA arms in a nanostar. This is done in the following way: First we compute the position, rcore, of the DNAns core. This is done by locating the base-pair (bp) closer to the FJC in each arm (three in total) and calculating the centre of mass (COM) using the position of the six nucleotides previously selected. Then, the unitary vectors (ê 1 ,ê 2 ,ê 3 ) pointing from the core of the molecule to the COM of the bp closer to the sticky end in each arm, define the direction of the three dsDNA arms.
To compute dp we find two vectors (A and B) that connect the end of any two arms, for example: These vectors define a plane whose normal can be found as the cross product n = A × B. Then, the distance dp from this plane From column two to five we show results at [NaCl]=0.15 M for n = 2, 3, 4 and 5, respectively. Top plots show the time evolution of the collective variable (dp). In each plot, two curves with different colour are shown. In turquoise we show results for a normal MD simulation without bias (γ = 1) and the other colour represents results from WTMetaD simulations with γ = 32. Bottom depicts the error analysis on dp as function of the block size. Insets show the temporal evolution of the Gaussian Height for the WTMetaD to the core of the molecule, is given by the projection of any of the three dsDNA arms (for exampleê 2 ) onto the unitary normal vector: The range of our CV is from dp = 0 (for a completely planar DNAns) to dp 1. The larger the value of dp the less planar the molecule. Note that dp cannot be equal to 1 due to the excluded volume interaction between arms.

The oxDNA model
We used the most recent implementation 1 of the oxDNA model into LAMMPS 2 (Large Scale Molecular Massively Parallel Simulator). Briefly, this model describes DNA at the single nucleotide level by means of a rigid body with additive-pairwise interaction sites. The potentials involved in the interactions accurately represent: the hydrogen bonding between complementary bases, the connectivity of the sugar-phosphate backbone, the excluded volume between nucleotides and also the stacking, coaxial-stacking and cross-stacking interactions. The relation between one simulation unit (SU) in the oxDNA and the international system (SI) units is: mass (M = 100AM U = 1.66 × 10 −25 kg), temperature (T = 3000K), length (Ls = 8.518×10 −10 m), energy ( s = k B T = 4.142 × 10 −20 J) and force (F = s/Ls = 4.863 × 10 −11 N ). The simulation time τ LJ = Ls M/ s = 1.7ps, comes naturally from the above quantities and it is employed to define a constant integration timestep δt = 0.001τ LJ .

Implementation
Here we use the plumed 3 library implemented in LAMMPS to perform Well Tempered Metadynamics of a three-armed DNA nanostar simulated with the oxDNA model. The sequence of this molecule is given in Table 1 of the main text and it corresponds to the one with n = 2 unpaired nucleotides at the FJC. We also performed WTMetaD to molecules with n = 3, 4 and 5 (see Fig. S1). The parameters used in the WTMetaD simulations are: γ = 32, B = 0.5k B T , ∆ = 0.05σ and τg = 4 × 10 5 timesteps. Each simulation was run for a total of 1.6 × 10 9 timesteps.
As described above, by choosing carefully the values of the bias factor (γ) in the WTMetaD we can compare results from different methods, namely, MD (γ = 1) and WTMetad (γ = 32). Top plots in Fig. S1 show the time evolution of dp. We observe that in simulations without the bias potential (shown in turquoise) dp does not explore the full range of possible values (dp is defined in the interval [0,1]). The maximum value reached is ∼ 0.65, it only happens briefly and for n = 5. Instead, in MetaD the CV is able to explore the whole range of configurations. Bottom plots show the error block analysis of the FEL, which is used to assess the convergence of the WTMetaD simulations. As expected, the error increases with the block size until it reaches a plateau in correspondence of a dimension of the block that exceeds the correlation between data points.

II Coarse-Grained Molecular Dynamics Simulations
Molecular dynamics simulations are performed using the model described in the main text. Here we refer to the details on those simulations. The Langevin integration of the system was carried out using LAMMPS in an NVT ensemble by a standard velocity-Verlet algorithm with integration time-step δt = 0.01. The position of particles in the system obeys the following equation where m = 1 is the mass of the particle, r represents its position. ξ is the friction and Λ(t) is the white noise term with zero mean which satisfies Λα(t)Λ β (s) = δ αβ δ(s − t) along each Cartesian coordinate represented by the Greek letters. The total potential field experience by a particle is U and depends on the interactions set for a particular system. In our coarse-grained model, each DNA nanostars is represented by a rigid body made up by ten particles. Seven beads account for the DNAns core and their excluded volume is modelled via a truncated and shifted Lennard-Jones (LJ) potential: if r < 2 1/6 σ, and U LJ (r) = 0 otherwise. Here σ = 2.5 nm represents the diameter of a spherical bead, = 1.0 parametrises the strength of the repulsion and r is the Euclidean distance between the beads. Patches are placed at a distance of 2.5σ from the core of the molecule and on the surface of the outermost bead along each arm. The sticky ends interaction, responsible for holding two DNAns together, is modelled by a Morse potential: The hard-sphere repulsion between beads forces patches to remain at a minimum distance of 0.12 σ. Right panel shows the plot of the morse potential used in simulations to set the attraction between patches (obtained from Eq. S7 with m = 25.0k B T , α 0 = 14σ −1 and r 0 = 0). Although the equilibrium distance between patches is set to zero, the minimum distance allowed by the excluded volume of the beads sets the energy Um(r = 0.12σ) ∼ 10 k B T (red dashed lines). for r < Rc. Here, r represents the distance between patches of two adjacent nanostars, r 0 = 0 is their equilibrium distance and Rc = 0.2σ is the cut-off distance of attraction. The amplitude of the potential is set to m = 25.0k B T and α 0 = 14σ −1 controls the width of the potential. These parameters were chosen to ensure that during the simulation any of the three arms of one ns can hybridize with any (but only one) of the arms of another ns. Importantly, the centre of a patch is placed at the edge of the last bead forming a dsDNA arm, i.e., at a distance of 0.5 σ from the centre of that bead. Since the hard-sphere repulsion of beads (set by the LJ potential in Eq. S6 ) covers a radius of (2 1/6 σ)/2 = 0.56σ, we expect hybridized patches to be at a minimum distance of ∼ 0.12σ. Therefore, the binding energy felt by hybridized patches is ∼ 10k B T (see Fig.S2), in agreement with the activation energy Ea mentioned in the main text. It is worth mentioning here that in the CG model we quote time in units of τ Br , the Brownian time, which is proportional to the time needed for a DNA bead to diffuse its own size. Finally, we note that the total excluded volume (V 1 ) of one DNAns is given by the sum of the volume of all its structural beads: In all our simulations the initial configuration is obtained by placing N nanostars in a cubic simulation box of length L = 40σ and with periodic boundary conditions. The volume fraction of this system is ρ = N V 1 /L 3 . We set m = 0 so nanostars cannot hybridized and we equilibrate the system for 5×10 5 τ Br . We then turn on the attraction between patches and record the evolution of our system for 10 6 τ Br . Results for planar and non-planar nanostars at different temperatures are shown in Fig. S3.

Melting curve
In the main text we show the melting profiles ϑ as a function of the temperature for the system at ρ = 0.01 (Figure 3(a)). This is obtained by averaging over the last 2 10 5 τ Br of the trajectories, the number of contacts Nc(t) (multiplied by the normalization factor 2/N f ). Examples are shown in Figs. S3(a,e). It is clear that all the samples start from an equilibrated configuration of unconnected nanostars (Nc(0) = 0); after turning on the attraction between patches the system evolves to a new steady state where a network is formed. The higher the temperature the more unconnected the stars and the smaller the value of < ϑ >.

The autocorrelation function
From the trajectory of Nc at long times, it is also possible to compute the time for network reconfiguration (τc). This is done by first finding the autocorrelation function of Nc(t), which in the discrete case has the form: where t f is the total number of time data points, t is the lag time and µ is the mean of the data. Results at different temperatures are shown in Figs. S3(b,f). The autocorrelation function follows an exponential decay c(t) ∝ e −t/τc , from which τc is computed.
In the main text we show τc as function of T in Figure 3(b).
We also note that, as expected, the data at lower temperatures has larger correlations. This is in agreement with the MSD plots shown in Figs. S3(c,g), where the system at lower temperatures has a smaller mobility.

IV Structural analysis
The radial distribution function The structure of our networks is analysed by computing the radial distribution function (RDF) of beads at the core of the DNA nanostars: where ρav = N/L 3 is the averaged density of core beads in the simulation, r ij represents the distance between the cores of molecules i and j. The sum counts the number of DNAns cores that are at a distance r.
In our simulations at fixed ρ = 0.01 there is a total of N = 175 nanostars (and core beads). We compute g(r) using the wrapped coordinates in the following way: (i) First, we choose the M core beads that are inside an sphere of radius R select = L/4 with centre at the origin of the box. (ii) We loop over each one of this M particles, and we count the number of core beads neighbours that are at a distance r. This search radius can cover up to a distance L/4 and in the count we consider all the core beads in the simulation (not only the M selected in (i)). (iii) We use this information in Eq. S10 to compute g(r). It is worth noting here that we repeat this procedure to compute g(r) in all the configurations in the steady state (at long times and when ϑ(t) has reached a plateau) and we take the average over these configurations. Figs. S3(d,h) show the average RDF obtained at different temperatures for planar and non-planar molecules, respectively. In both cases, g(r) shows a maximum located near to 5σ, corre- DNAns) as function of the temperature for planar (dp = 0, in blue) and non-planar (dp = 0.6, in orange) molecules. All the data shown here correspond to simulations performed at ρ = 0.01 sponding to the average distance between the cores of two bound nanostars. As temperature increases, the hight of this maximum decreases, meaning that there are less contacts between nanostars. Figures S4 shows the size of the largest cluster in the networks (normalized by the number of nanostars N=175) made of either, planar or non-planar molecules. Over the whole range of temperatures explored, the size of the largest cluster (Max(Cs)) is larger in the case of networks made out of planar molecules.

V Varying d p
In the main text we show the fraction of connected DNAns ϑ and the relaxation time τc as function of the planarity (dp), when simulations are performed at a fixed temperature (  and concentration (ρ = 0.01). Here we provide the full trajectories from which these results were obtained. Fig. S5(a-b) show the temporal evolution of Nc(t) and c(t). We observe that as dp increases, more contacts between DNAns are formed and the contacts-correlation is larger. The RDF in Fig. S5(c) shows that the two local maxima located at r ∈ [6, 8] are developed when dp reaches values ∈ [0.5, 0.6]. At higher values of dp, the amplitude of g(r) at these local maxima shrinks, and the amplitude of the global maximum (at r close to 5) increases. This is consistent with the increase in the number of box-like structures observed in the network diagrams of Fig. 5 in the main text. In the extreme case dp = 0.9, one local maxima appears at particular small values of r < 5. This is in agreement with the structure made by four DNAns and depicted in the snapshot of Fig. 5(d) in the main text and for which DNAns connect with a bending angle φ ∼ 20°, as shown in Fig. 6(c).

VI Green-Kubo relations
The autocorrelation function (G(t)) is defined as: where P αβ represents the out-off diagonal component (Pxy, Pxz and Pyz) of the stress tensor. Results from long simulations (1.6 × 10 6 τ Br ) are shown in Fig. 4(c) of the main text. The autocorrelation was computed using the multiple-tau correlator method described in reference 4 and implemented in LAMMPS with the fix ave/correlate/long command. This computation ensures that the systematic error of the multiple-tau correlator is always below the level of the statistical error of a typical simulation (see LAMMPS documentation). The viscosity (η) of the system is then obtained by computing the following integral in practice we compute the numerical discrete integral of the curves in Fig. 4(c). We obtained the viscosity in simulation units for planar (ηp = 573 k B T τ Br /σ 3 ) and non-planar (ηnp = 29 k B T τ Br /σ 3 ) networks at ρ = 0.01. To convert these results to real units we use σ = 2.5 nm, k B T = 4.11 pN · nm and the Brownian time: (S13) Assuming that the solvent is water (with viscosity ηs = 1 mP a · s) we obtain τ Br = 36 ns. Therefore, the viscosities of the two networks map to ηp = 5.5P a · s and ηnp = 0.28P a · s.

VII Varying concentration
Up to now, we have presented results from simulations performed at a constant volume fraction ρ = 0.01. In this section, we investigate the effect of increasing concentration for networks made of planar and non-planar nanostars. Without loss of generality we choose dp = 0.6 for the non-planar case. First, we compute the size of the largest connected component of the network (Max(Cs)). When Max(Cs)= N , there is only one cluster in the system that is formed by the connection of all DNAns. In contrast, when the value of Max(Cs)< N , nanostars are split into clusters. In Fig. S6(a) we report the value of Max(Cs) obtained from simulations at different volume fractions, note that the results are normalized by the number of nanostars at each concentration. We observe that as ρ increases, the fraction of nanostars connected to the largest cluster also grows, for both planar and non-planar molecules. However, the value of Max(Cs) is consistently lower in the non-planar case. At ρ = 0.06, most of the nanostars are part of a single cluster regardless their planarity. In the following, we focus in results obtained at this value of ρ, but the same general results were obtained at different volume fractions.
In Fig. S6(b) we corroborate that at a fixed volume fraction (either ρ = 0.01 or ρ = 0.06), the equilibrium value of the fraction of contacts (θ) is larger for non-planar nanostars. The RDF is shown in Fig. S6(c), we observe that although the height of the peaks has decreased with the concentration, the region at which the maxima appear is conserved (compare blue and purple, or orange and red). Finally, Fig. S6(d) shows the network diagrams obtained for both designs. For the non-planar case a few small clusters appear, indicating a slightly lower connectivity than the planar case. The histogram depicting the frequency of the nanostars with certain degree of connection is also shown. As it can be seen, the number of DNAns fully connected (deg(v) = 3) is larger for the non-planar molecules. This is all in agreement with the discussion in the main text, where we show results mainly for ρ = 0.01.

VIII Non-rigid model
In the discussion so far, it is clear that the planarity of the DNAns plays a major role in determining the geometry and topology of DNA hydrogels, and in consequence, the planarity also determines their viscoelastic properties. However, one may wonder if this result could be somehow related to the fact that we use rigid bodies to represent DNAns in our coarse-grained model. Here we first introduce a coarse-grained model of DNAns with additional degrees of freedom, and then we use it to prove that our results hold when we allow small fluctuations of the angle α ij between the dsDNA arms of individual nanostars. The non-rigid model of DNAns is depicted in Fig. 7 and it can be seen as a generalization of the model used in the main text. It consists of four beads (depicted in red) forming the core of the molecule, and three beads (depicted in blue) with a patch attached to their edge (depicted in cyan). Each of these patchybeads is integrated as a rigid body. All beads have an excluded volume of 1σ modelled with the LJ potential in Eq. S6. The attraction between patches is modelled with the Morse potential in Eq. S7. Neighbouring beads are connected via a harmonic potential: where r 0 = 1.1224σ is the equilibrium distance between beads and bond = 300k B T holds beads in place and only allows for small fluctuations around r 0 .
We use an improper dihedral (called umbrella style in LAMMPS) to have control over the planarity of the DNAns. The potential representing this interaction has the form: where w is the angle between the IL axis and the IJK plane in Fig. 7(b). imp = 300k B T is the stiffness of the potential imposing an equilibrium angle w 0 . For planar molecules we set w 0 = 0°and for non-planar molecules we set w 0 = 87°, which corresponds to dp ∼ 0.6. The angle interactions are sketched in Fig. 7(c): α represents the angle between any two arms of a nanostar, β is the angle used to align the beads forming a single dsDNA arm, and γ is the angle setting the alignment of the patch with the dsDNA arm. Note that a total of nine angles are needed (3 of each type) to represent a single DNAns. All the angle interactions are modelled by harmonic potentials with equilibrium values α 0 = 120°, β 0 = 180°a nd γ 0 = 180°to resemble the structure of a DNAns. Because each of the dsDNA arms is 20 − 30 base-pairs long, smaller than lp = 150 bp (the persistence length of DNA), we expect arms to be straight, and therefore, we set β = 300k B T . We decided to set γ = 300k B T to allow for small fluctuations of the alignment of the patches with the arm. Recall that this choice provides a natural distribution of bending angles (φ) at the place of hybridization between two DNAns (see Fig. 6 of the main text). To study the impact of the fluctuation of the angles between arms within a single DNAns, we use two different values of α. At high values ( α = 300k B T ), we expect a similar behaviour to the rigid analogues of nanostars studied in the main text. On the other hand, for small values of α, nanostars fluctuate around the Y-shaped conformation in a qualitatively compatible way with the conformations found in simulations with the oxDNA model for the planar design (with n = 2 at [NaCl]=0.15 M, see for instance Supplementary Movie 2). In practice, we performed simulations with the non-rigid model using w 0 = 0°(planar), and we decreased the value of α until we observed that for α = 6k B T , nanostars are able to deform into other conformations.
In Fig. S8 we show results comparing the different models: rigid (used in the main text), non-rigid with α = 300k B T and non-rigid with α = 6k B T . As we expect, the two first models produce outcomes that are very similar. This is particularly truth for the planar DNAns, where the red and orange curves are very close in Panels (a)-(d). For the non-planar case, curves in blue and cyan deviate a little bit. We attribute this to the fact that w 0 = 87°doesn't provide exactly the value of dp = 0.6. We also note that for the non-rigid model the RDF shows a shift of the position of the local maxima.
When we compare the results for the non-rigid model with the same planarity, the case with α = 6k B T is more dynamic. Unsurprisingly, the MSD shows that the mobility is lower when α = 300k B T compared to the result for α = 6k B T . This is also reflected in the autocorrelation of the stress tensor, where G(t) is lower for the latter case.
Finally, Panel (e) displays results for the non-rigid model with α = 6k B T and it shows that for non-planar molecules, the mobil- ity is higher (top) and the viscosity is lower (bottom), than those of the networks formed by planar molecules. This supports the results obtained with the rigid model and allows us to be confident that allowing fluctuations of the angle between arms in a single DNAns, will not change the general results found in the main text. It would also be interesting to perform additional simulations (with higher resolution models of DNA) using the methods described in this manuscript, in order to tune the rest of the potentials in our non-rigid model.