Mechanistic insight into the structure, thermodynamics and dynamics of equilibrium gels of multi-armed DNA nanostars

Abstract

The unique sequence specificity rule of DNA makes it an ideal molecular building block for constructing periodic arrays and devices with nanoscale accuracy and precision. Here, we present the self-assembly of DNA nanostars having three, four and five arms into a gel phase using a simplistic coarse-grained bead-spring model developed by Z. Xing, C. Ness, D. Frenkel and E. Eiser (Macromolecules, 2019, 52, 504–512). Our simulations show that the DNA nanostars form a thermodynamically stable fully bonded gel phase from an unstructured liquid phase with the lowering of temperature. We characterize the phase transition by calculating several structural features such as the radial distribution function and structure factor. The thermodynamics of gelation is quantified by the potential energy and translational pair-entropy of the system. The phase transition from an arrested gel phase to an unstructured liquid phase has been modelled using a two-state theoretical model. We find that this transition is enthalpy driven, and loss of configuration and translational entropy is counterpoised by enthalpic interaction of the DNA sticky-ends, which gives rise to a gel phase at low temperature. The absolute rotational and translational entropy of the systems, measured using a two-phase thermodynamic model, also substantiates the gel transition. The slowing down of the dynamics upon approaching the transition temperature from a high temperature demonstrates the phase transition to a gel phase. A detailed numerical simulation study of the morphology, dynamics and thermodynamics of DNA gelation can provide guidance for future experiments, is easily extensible to other polymeric systems, and is expected to help in understanding the physics of self-assembly.

---

I. Introduction

Utilizing the unique base-pairing specificity rule, DNA can be exploited as an ideal building block to design highly organized supramolecular materials with nanoscale accuracy.¹ Purposely sequenced single-strand DNA bases can form DNA motifs via reciprocal exchange between adjacent duplexes. These motifs, via hierarchical self-assembly, form a wide range of highly organized artificial nanostructures. The concept of DNA nanotechnology has its root back in the 1980s when Ned Seeman of New York University first designed an immobile four-way DNA junction using four different single-stranded DNA.^2,3 This idea led to the foundation of the structural DNA nanotechnology field to construct synthetic two-dimensional and even three-dimensional DNA nanostructures.^4,5 Later, the “origami” method revolutionized the field by enabling the creation of more complex and larger DNA nanostructures.⁶ Since then, the field has expanded rapidly due to its potential applications in diverse fields such as nanomedicine, nanoelectronics, synthetic biology, and nanorobotics.^7–13

DNA can be assembled into complex 1D, 2D and 3D DNA nanostructures with precise arrangements of helices either via DNA origami or DNA tile techniques. In a similar fashion, small DNA nanostructures can self-assemble into a complex polymeric network at desired temperature forming hydrogels.^14–18 The resulting gels formed via hydrogen bonding and other dispersive interactions can be classified as physical hydrogels and are biodegradable, bio-compatible, and non-toxic. They can encapsulate drugs, cargo, and even cells with extensive applications such as controlled drug delivery, 3D cell culture, bio-printing, cell transplant therapy, tissue engineering and other biomedical applications.^14–17,19 The mechanism of DNA hydrogel formation is relatively simple. Single-stranded DNA first hybridize to form small DNA nanostars with sticky ends and then these nanostars exploiting the sticky-end cohesion bind with each other to form a spanning polymeric network. By tuning the length and sequence of sticky ends, it is even possible to control the gel formation thermally. Nagahara and Matsuda first reported a hybrid polyacrylamide–DNA hydrogel.²⁰ Later Liu, Zhou, and co-workers using a one-pot step-wise self-assembly protocol, prepared the first pure DNA supramolecular hydrogel.²¹ In 2006, Luo et al. fabricated a hydrogel under physiological conditions entirely from branched DNA by enzymatic ligation, which exhibited amazing shape memory properties.¹⁴ Several other experimental observations have been reported on different aspects and applications of DNA hydrogels.^19,22–32 More recently, DNA based hydrogels have been utilized to stimulate membrane endocytosis, which leads to enhanced cell spreading and invasion of cells.¹⁹ Notably, Eiser et al. extensively probed the microrheological and viscoelastic properties of DNA hydrogels.³³ The collective phase behaviour and dynamics of DNA hydrogels have been investigated both experimentally and numerically.^30,34–44

One of the critical issues of the DNA hydrogel paradigm is predicting the complex three-dimensional self-assembled structures. The kinetics of the multidimensional aggregation process are relevant to solve many fundamental physics issues. Also, how the morphology of constituent DNA nanostars affects the disordered arrested gel states is known poorly. Most of the recent studies focus mainly on the synthesis, modelling and applications of DNA hydrogels and the fundamental physics of gelation still lacks a good understanding. Therefore, the challenge to microscopically understand the underlying physics of self-assembly remains.

All-atom models, which were found to be very effective in addressing self-assembly and mechanical properties of DNA nanostructures, are too expensive to study DNA hydrogel formations.^45–48 Coarse-grained (CG) models with different levels of descriptions have found to be more effective in exploring different aspects of gelation.^{35–39,43,49} Although, most of the base-pair level coarse-grained models are too detailed and computationally demanding and are thus often not useful to study the self-assembly of DNA nanostars.^50–55 Recently, a simple bead-spring model has been introduced by Xing et al. in which the sticky ends are replaced with a single patch.⁵⁶ The model has captured the self-assembly and microrheological properties of gelation of three armed “Y” shaped DNA nanostars.

In this work, using the simple bead-spring model we have studied the structure and thermodynamics of gelation of DNA nanostars having more complex morphologies. We have modelled three different classes of DNA nanostars: (i) trivalent nanostars (referred to as Y-DNA), (ii) tetravalent nanostars (referred to as X-DNA) and (iii) pentavalent nanostars (referred to as 5WJ). In the methodology section, a brief description of models and methods implemented in this study is given. In the result section, we have first discussed the structure of DNA hydrogel, followed by thermodynamics of gelation, where we have explicitly evaluated the potential energy, translational pair entropy and change in enthalpy and entropy upon gelation. Then, we employed the two-phase thermodynamic method (2PT) to calculate the absolute translational and orientational entropy to distinguish different phases of hydrogel. Finally, the dynamics are probed by velocity autocorrelation and mean square displacement of the systems.

II. Models and numerical methods

A. Coarse-grained model

We employed a bead-spring CG model to represent the structure of DNA nanostars which was originally developed by Eiser and co-workers.⁵⁶ All the force field (FF) parameters like mass (m_LJ), length (σ_LJ), energy (ε_LJ), and the Boltzmann constant (k_B) are expressed in reduced Lennard-Jones (LJ) units. In terms of these quantities, the reduced time is defined as . The CG structure is composed of two types of charge-less beads–one large bead to replicate the double-stranded nature of DNA and one small patch to represent the sticky ends (Fig. 1). Two types of small patches are defined (designated as patch-I and patch-II). An attractive potential is implemented between the opposite patches only, while a repulsive potential is applied between the same type of patches. This is implemented to prevent any three-body interaction between the patches. The non-bonded excluded volume interaction between the beads is modelled using truncated and shifted LJ potential, known as Weeks–Chandler–Andersen (WCA) potential truncated at the minima of the potential (distance 1.12σ_LJ).where r is the distance between two particles, ε is the depth of the potential and σ is the size of the particle. U′ is set in such a way that U_WCA(r = 1.12σ_LJ) = 0. In our simulation, we set ε = ε_LJ, σ = σ_LJ and r_cutoff = 1.12σ_LJ for the interaction between two large beads. This gives rise to a repulsive potential which prevents any unwanted overlap between the arms of the nanostars. For non-complementary patches, we set ε = ε_LJ, σ = 0.67σ_LJ and r_cutoff = 0.67σ_LJ which leads to a repulsive interaction between the same patches and prevent multiple attraction between different patchy ends. Suppose a bond between patch-I and patch-II is formed, then the repulsive potential prevents any possible bond formation of patch-I–patch-I or patch-II–patch-II which gives each arm a strict valency of 1 (Fig. 1(d)). The attractive interaction between the complementary patches is described using full LJ potential truncated at a longer distance as given using eqn (2),For patch-I and patch-II interaction, we set ε = 4.5ε_LJ and σ = 0.2σ_LJ. This leads to an attractive potential at a very short distance (Fig. 1(e)), which mimics the sticky-end attraction of different DNA arms.

Fig. 1 Coarse-grained models of the DNA nanostars investigated in this study. (a) Y-shaped DNA nanostars, (b) X-shaped DNA nanostars and (c) 5 armed (5WJ) DNA nanostars. (d) Non-bonded interactions between different CG beads. (e) Pairwise LJ potential between different beads. (f–h) Initial configuration of the system with different densities. The dimension of each box is 30 × 30 × 30 σ_LJ³. The number of nanostars in each simulation box is (f) N = 50, (g) N = 150, and (h) N = 250, respectively.

All the neighbouring beads are connected by a harmonic potential of the form,

U_bond = K_r(r − r₀)² (3)

where K_r is the spring constant and r₀ is the equilibrium bond length. The value of K_r is set to a very high value of 300ε_LJ/σ_LJ² to allow only small fluctuation in the bond length and maintain the structural integrity of the nanostars. The value of r₀ between two large beads is set to 0.96σ_LJ, and between the large bead and any patchy particle is set to 0.56σ_LJ.

Similarly, the angle bending potential involving three beads is given using the following harmonic potential,

U_angle = K_θ(θ − θ₀)² (4)

where K_θ is the spring constant and θ₀ is the equilibrium bond angle. Like K_r, the value of K_θ is also set to a very high value of 300ε_LJ/radian² to maintain the Y-shape geometry. The equilibrium bond angles involving the central bead are θ₀ = 120° for Y-shaped nanostars, 90° for X-shaped nanostars, and 72° for 5WJ nanostars. The equilibrium bond angle involving three beads in the branch of the nanostar is 180°.

B. System and simulation details

The schematic illustration of the coarse-grained model of the DNA nanostars (Y-shaped, X-shaped and 5WJ) is shown in Fig. 1. For each DNA nanostar, we built three systems with 50, 150 and 250 units. We placed the nanostars in a periodic box with random positions and orientations. The dimensions of the periodic box were 30 × 30 × 30 σ_LJ³. The density, volume fraction and other details of the simulated systems are given in Table 1. The volume fraction is calculated based on the volume of each nanostar. The volume of the each arm is around ∼πr²l with r = 0.56σ_LJ and l = 2.48σ_LJ. Based on this calculation, the volume of one Y-shaped, X-shaped and 5WJ nanostar is 7.33σ_LJ³, 9.77σ_LJ³, 12.22σ_LJ³ respectively.

Table 1 Details of the studied systems

System	Number of nanostars	Number density (×10^−3)	Volume fraction (%)
Y-DNA	50	1.85	1.35
	150	5.55	4.05
	250	9.25	6.75
X-DNA	50	1.85	1.81
	150	5.55	5.43
	250	9.25	9.05
5WJ	50	1.85	2.26
	150	5.55	6.78
	250	9.25	11.30

---

We performed CG molecular dynamics simulation using the LAMMPS simulation package. The time evolution of the position of each bead was performed by employing the Langevin equation of motion given as,

where r and m are the position and mass of a bead respectively. U_tot is the total potential energy given as U_tot = U_WCA + U_bond + U_angle. γ is the damping factor. The value of γ is set to a high value of 100 to mimic the over-damped condition of the gel phase. ξ(t) is the stochastic noise coming from the interaction of the particle with the heat bath, which can be written as, . R(t) is the delta-correlated Gaussian white noise. For each system, we varied the temperature from 0.075ε_LJ/k_B to 0.775ε_LJ/k_B with an interval of 0.025ε_LJ/k_B. For more details of the model and simulation details, readers are referred to the original paper by Eiser et al.⁵⁶ The time step of integration is set to 0.005τ_LJ. For each temperature, we run a total of 10⁷ steps. To ensure equilibration, we make sure that the energy and the number of bonds formed attain a stationary value.

III. Results and discussion

A. Structure

In order to see the formation of a gel phase, we have explored a wide range of temperatures and densities. At high temperatures (T* > 0.5), the nanostars are randomly oriented in a gas-like structure. As the temperature is decreased, the nanostars self-assemble into a complex percolating network (Fig. 2). In most cases, such a network percolates over the whole simulation box. To characterize the structure of such complex networks, we have calculated the radial distribution function (RDF),where N is the total number of nanostars, ρ is the average number density of systems, and 〈…〉 denotes the ensemble average. The sum counts the total number of pairs at a distance r. The RDF of the central bead of each nanostar is calculated using VMD software⁵⁷ and has been shown in Fig. 3. At high temperatures (T* > 0.5), the nanostars form a gas-like unstructured fluid and the RDF is almost flat with no peaks. With the lowering of temperature (T* < 0.3), the nanostars form networks and we see the emergence of clear peaks in the RDF. The position of the first peak in the RDF can be estimated from the geometric criteria of the bonded nanostars. The length of an arm of a nanostar is around 2.48σ_LJ. When two nanostars are bonded, the distance between their central beads is around 4.96σ_LJ. The location of the first peak of the RDF is also at a separation of 4.96σ_LJ. Similarly, the location of the 2nd and 3rd peaks can also be described by using different arrangements of nanostars in the network (Fig. 3(d), (h) and (l)). Also with increasing densities, we observe a decrease in the value of the peak height of the first peak of the RDF. As the density increases the number of available pairs between r to r + dr remains almost the same but the total number of nanostars, N_pairs and number density, ρ increase. As a result, the height of the first peak of g(r) decreases.

Fig. 2 Structure of a DNA nanostar mixture at different temperatures. Y-DNA at temperature (a) T = 0.775ε_LJ/k_B and (b) T = 0.075ε_LJ/k_B. X-DNA at temperature (c) T = 0.775ε_LJ/k_B and (d) T = 0.075ε_LJ/k_B. 5WJ DNA at temperature (e) T = 0.775ε_LJ/k_B and (f) T = 0.075ε_LJ/k_B. The number of nanostars in each system is, N = 50.

Fig. 3 Radial distribution function (RDF) of the nanostar mixtures at different temperatures and densities. (a–d) RDF of Y-shaped nanostar mixtures. The number of nanostars in the systems is (a) N = 50, (b) N = 150, and (c) N = 250. (d) Location of the 1st and 2nd peak in the RDF when two Y-shaped nanostars are bonded. (e–h) RDF of X-shaped nanostar mixtures. The number of nanostars in the systems is (e) N = 50, (f) N = 150, and (g) N = 250. (h) Location of the 1st and 2nd peak in the RDF when two X-shaped nanostars are bonded. (i–l) RDF of 5WJ nanostar mixtures. The number of nanostars in the systems is (i) N = 50, (j) N = 150, and (k) N = 250. (l) Location of the 1st and 2nd peak in the RDF when two 5WJ nanostars are bonded.

Another interesting fact is that the network structure is amorphous.³⁰ However, from the RDF calculation, it is not clear whether the structures show crystalline or noncrystalline phases. To justify the absence of crystalline ordering in the network, we have computed the structure factor (SF), SF(q). The SF of central beads has been calculated by taking the Fourier transform of the RDF,

where q is the wave vector, r_i is the position of the particle, and 〈⋯〉 denotes the ensemble average. The structure factors of the central bead of various nanostar mixtures are shown in Fig. 4. Similar to the RDF, the SF(q) shows a flat profile indicating that mixtures form an unstructured fluid at high temperatures (T* > 0.5). At low temperatures (T* < 0.3), DNA nanostars undergo structural ordering. However, we do not see any signs of crystallization. The absence of any Bragg's peak in the SF(q) substantiates the amorphous nature of the network. Indeed, in several previous experimental and theoretical studies, it has been reported that the gels of DNA never crystallize.^30,43 Our calculation of the SF(q) reaffirms the observation. At very low T, we observe a trihedral network for Y-DNA mixtures (observed from the number of peaks in SF(q)), whereas for X-DNA we see a tetrahedral order. For 5WJ DNA nanostars, we find a pentagonal order at very low T. At high densities, i.e., when the number of nanostars in the system is N = 250, the trihedral/tetrahedral/pentahedral network of the different nanostar systems becomes more prominent. Also, at lower densities, we observe an increase in SF(q) as q → 0. This is an indication of high compressibility at low density. At low density, liquid-like characteristics are more dominant. As a result, we observe high compressibility. However, as we increase the nanostar density, gel characteristics start to dominate and the compressibility reduces. Similar observations were reported by Sciortino et al. in their study of tetra-valent DNA hydrogels.⁴³

Fig. 4 Structure factor (SF) of the nanostar mixtures at different temperatures and densities. (a–c) SF of Y-shaped nanostar mixtures. The number of nanostars in the systems is (a) N = 50, (b) N = 150, and (c) N = 250. (d–f) SF of X-shaped nanostar mixtures. The number of nanostars in the systems is (d) N = 50, (e) N = 150, and (f) N = 250. (g–i) SF of 5WJ nanostar mixtures. The number of nanostars in the systems is (g) N = 50, (h) N = 150, and (i) N = 250.

B. Thermodynamics

The formation of a network is mainly driven by the interaction between the patchy beads resulting in the formation of non-covalent bonds, which resemble the hydrogen bonds between the sticks-ends. At low temperatures, this non-covalent bond formation in the gel phase gives rise to a large gain in enthalpy. To have a quantitative understanding, we have computed the potential energy (PE) of a nanostar mixture at different temperatures. The potential pair energy is calculated by summing Lennard-Jones potential energies between all the non-bonded atom pairs. In Fig. 5, we have plotted the PE per atom for all the nanostars at different temperatures. We find that the PE is positive when the temperature is high. PE becomes negative at low temperatures. Also, as we increase the temperature we observe a discontinuity in the slope of PE over a small range of temperature (around T* ∼ 0.3 to 0.5). The nanostars undergo from a fully bonded state to an unstructured state within this temperature range. The discontinuity also implies that the systems go through a continuous phase transition as the temperature is increased or decreased. This morphological transition from a fully bonded state to an unstructured fluid can be well described using a two-state theoretical model as discussed below.

Fig. 5 Potential energy (PE) per atom of the nanostar mixtures at different temperatures and densities. The systems consist of (a) Y-shaped, (b) X-shaped and (c) 5WJ nanostars. The fraction of intact bonds at different temperatures and densities. The systems consist of (d) Y-shaped, (e) X-shaped and (f) 5WJ nanostars. (g) Change in entropy and (h) change in enthalpy as the system goes from a gel state to an unstructured state.

To define the non-covalent bond formation between two nanostars or association between two nanostars, we define a geometric distance cutoff of 0.5σ_LJ. Since the non-bonded PE between the patches is negative when they are within this confining cutoff distance, our definition of a bonded network is reasonable. We calculate the total number of bonds in each system and divide it by the maximum number of bonds formed at very low temperatures to get the fraction of connected bonds. Similar to the behaviour of PE, we observe that the system goes from a state where none of the nanostars is bonded to a state where all the nanostars are bonded and form a percolating cluster (Fig. 5(d)–(f)). Several experimental and numerical results on thermo-reversible gel transition show a similar behaviour.^{30,34–38,43} Also, we find that this behaviour is almost independent of the densities of the system. This sharp phase transition has a significant contribution of free energy for nanostar association, which has primarily two parts–entropic and energetic contribution. We try to estimate the free energy change due to the formation of bonds using a simple two-state model. According to the model, the fraction of bonds, f, can be written as,

where ΔE is the change in enthalpy and ΔS is the change in entropy upon the phase transition. We fitted eqn (8) for the fraction of bonds with ΔE and ΔS as the fitting parameter (Fig. 5 (d)–(f)). We find that the entropic change is exactly the same as the enthalpic change upon gel to unstructured fluid phase transition (Fig. 5(g), (h) and Table 2). The entropic change originates from the localization of the patchy sites in a small confining region when bonded compared to the unbonded state where nanostars have full orientation and translation freedom. The contribution of enthalpy change comes from the bond formation between the patchy beads in the gel phase. This enthalpy change due to the formation of bonds is comparable to the depth of attractive LJ potential of the patchy beads which in our model is equivalent to 4.5ε. The entropic loss in the gel phase is compensated by the enthalpic gain. This large enthalpic gain due to patchy bead interaction is primarily responsible for this phase transition from unstructured liquid to a gel phase. Next, we try to estimate how much translational entropy contributes to this phase transition. The total pair distribution function (g(r,θ)) of a system can be decomposed into two parts– the radial pair distribution function (g(r)) and orientational pair distribution function (g(θ|r)),

g(r,θ) = g(r)g(θ|r) (9)

Table 2 Entropy change (k_B) between unstructured fluid and gel phases, calculated using different methods^a

System	Number of nanostars	Gel-to-sol transition temperature	Enthalpy change from two-state model (εLJ)	Entropy change from two-state model (kB)	Translational pair entropy change (kB)	Absolute translational entropy change (kB)	Absolute rotational entropy change (kB)
a Entropy difference between the gel phase and fluid phase for the three former methods is evaluated by taking the difference in average entropy between the two states. Average entropy is calculated by averaging extreme five state points of entropy vs. temperature graphs.
Y-DNA	50	0.40	4.36	10.84	0.92	5.95	3.89
	150	0.41	3.92	9.51	0.38	5.78	3.64
	250	0.43	3.72	8.69	0.25	5.70	3.62
X-DNA	50	0.40	5.11	12.73	3.16	10.94	5.08
	150	0.40	3.49	8.72	1.07	9.64	4.97
	250	0.41	3.12	7.57	0.73	9.42	4.97
5WJ	50	0.43	4.92	11.55	2.72	10.02	4.85
	150	0.44	4.01	9.20	1.03	9.23	4.92
	250	0.45	3.71	8.21	0.74	9.01	4.97

---

For a complex molecule like DNA nanostars, the orientational pair distribution function and consequently the orientational pair entropy is difficult to estimate. However, the translational pair entropy can be estimated from the radial distribution function using the following relation,⁵⁸

The calculated translational pair entropy at different temperatures shows a similar behaviour as potential energy (Fig. 6). When the nanostars are bonded and form a gel phase, they lose many degrees of freedom and the translational pair entropy of the system becomes significantly low (Table 2). However, as the temperature increases and the fraction of connected bond reduces, the system gains more degrees of freedom and translational pair entropy increases (Table 2). Also, we observe that with increasing density, the translational pair entropy for each nanostar system decreases. This is due to the fact that when the density of a system increases the available free volume fraction decreases. The reduction of the free volume fraction in turn reduces the translation motion of each nanostar. Therefore, the translational pair entropy of the system decreases with increasing density. The qualitative characteristic of the translational pair entropy change with temperature is similar for all the systems.

Fig. 6 Translation pair entropy of different nanostar systems at different temperatures and densities. The systems consist of (a) Y-shaped, (b) X-shaped, and (c) 5WJ nanostars.

It is now apparent that the system losses entropy as it transitions to the gel phase. In order to evaluate the absolute entropy of the system, we have employed the robust two-phase thermodynamic (2PT) model developed by Lin et al.^59–61 The 2PT method is built upon the simple hypothesis that the density of states (DOS) can be decomposed into solid-like (oscillatory) and gas-like (diffusive) components. First, for any poly-atomic molecule, the total DOS can be decomposed into the translational, rotational and vibrational parts,

S^tot_DOS(ν) = S^trn_DOS(ν) + S^rot_DOS(ν) + S^vib_DOS(ν) (11)

The DOS is computed from the Fourier transform of the velocity autocorrelation function (VACF). The density of state function S_DOS(ν) is defined as the mass-weighted sum of the atomic spectral densities. This can be obtained from the Fourier transform of VACF obtained from the simulation as follows,where N_atoms is the total number of atoms in the system, m_l is the mass of the l-th atom and v^k_l is the velocity of the l-th atom in the direction k(x,y,z). From the translational component of the center of mass velocity, we calculate the S^trn_DOS. Similarly, from the angular and vibrational velocity, we calculate the S^rot_DOS and S^vib_DOS, respectively. The total DOS of each nanostar system in different frequency regimes is plotted in Fig. 7(a) and (b). For all the systems, we find that the DOS has a pronounced peak in a low-frequency regime whose height increases with increasing temperature. At high temperatures, the DOS exponentially decays, which indicates that the structures are more gas-like. However, at low temperatures, the DOS first increases and then decreases. This DOS in a low temperature regime is similar to the DOS of liquid systems, which indicates that some structural ordering has emerged in the system. Also, in a high-frequency range, the DOS has several peaks. These peaks correspond to different translational, orientational stretching and vibrational modes of the nanostars. To validate that the DOS correctly captures the dynamics of the system, we computed the diffusion constant of the system from the zero frequency DOS using the following relation, , where D is the diffusion constant.^59,60 We observe that the diffusion constant values are very close to the diffusion constant obtained from the slope of mean-squared displacement versus time data. More details about the diffusion constant calculation are in the next section.

Fig. 7 Density of states (DOS) of different nanostar systems at three different temperatures. (a) Total DOS for all frequency regimes and (b) DOS in the low-frequency regimes.

In the 2PT model, each part of the DOS is decomposed into solid-like and gas-like contributions.

S^α_DOS(ν) = Ψ^α_solid(ν) + Ψ^α_gas(ν); α = trn/rot/vib (13)

Then the entropy (S) or any thermodynamic quantity is determined from the DOS by introducing the weight factor W_solid and W_gas,Using the 2PT method, we next evaluated the absolute translational entropy (S_{translational}) and rotational entropy (S_rotational) for different nanostar systems (Fig. 8). First, the LJ units are converted into real units using Argon interaction parameters and used in the 2PT code. Once entropy values are obtained using the 2PT code, results are again mapped to the LJ units. The computed S_{translational} and S_rotational show a similar behaviour to the other thermodynamics parameters like PE and transitional pair entropy. As the temperature increases, both S_{translational} and S_rotational increase with temperature and at around temperature 0.4ε_LJ/k_B, there is an abrupt jump in the entropy. Also, the systems lose more translational entropy compared to the rotational entropy in the gel phase (Table 2). Another interesting fact about the gel system is that unlike liquids they lose huge rotational entropy as the temperature is decreased. And this specific property makes it distinguishable from the conventional liquid phase.

Fig. 8 Absolute translational and rotational pair entropy per nanostar for different systems at different temperatures and densities. Translational entropy of the (a) Y-DNA, (b) X-DNA, and (c) 5WJ system. Rotational entropy of the (d) Y-DNA, (e) X-DNA, and (f) 5WJ system.

C. Dynamics

To better understand the structure and thermodynamics of the gelation with the dynamical arrest of the system, we next computed the mean square displacement (MSD) as a function of time for different temperatures and densities. The MSD of the central bead of each nanostar is evaluated using the following relation,where N is the total number of beads, t is the time difference, t′ is the time origin and the angular bracket denotes an average over all the time origins. The calculated MSD of different systems is plotted in Fig. 9. The diffusion constant is calculated from the linear region of the MSD using the following relation,

MSD = 〈r²〉 = 6Dt (16)

where D is the diffusion constant of the system. At high temperatures, the MSD of the system increases linearly with time. Also, as the temperature decreases and the nanostars make bonds with one another, and the diffusion constant also decreases (Fig. 9). Moreover, when the system forms a complete percolating network and goes to a gel phase, the Einstein relation of the diffusion constant, i.e., eqn (16) no longer holds. Therefore, we are restricted to calculate the diffusion constant up to temperature, T* = 0.4. Also, we have calculated the diffusion constant from the density of states at zero frequency from the 2PT method,^59,60where C(t) is the mass-weighted sum of the velocity autocorrelation function. The diffusion constant calculated using the 2PT method is slightly higher compared to the diffusion constant calculated from MSD. As in MSD calculation, we only considered the central bead which removes any kind of rotational motion. However, in the 2PT method, the translational, rotational and vibrational motion of each nanostar is considered while calculating the velocity autocorrelation function. Thus, the diffusion constant in MSD is slightly underestimated. Also, we find that with increasing density, the value of D decreases for every system. This is expected as with increasing density, the accessible phase space volume of the nanostars decreases. Another interesting fact is that the diffusion constant for any particular density is the highest for Y-DNA and the lowest for 5WJ nanostars while the diffusion constant of X-DNA lies in between them. This can also be explained from the volume fraction of each system. The volume fraction of 5WJ is the highest compared to the other two nanostar systems, reducing the available phase space volume to diffuse and this is reflected in the overall diffusion of the system.

Fig. 9 Mean squared displacement (MSD) of different nanostar systems at different temperatures. (a–c) MSD of Y-shaped nanostar mixtures. The number of nanostars in the systems is (a) N = 50, (b) N = 150, and (c) N = 250. (d–f) MSD of X-shaped nanostar mixtures. The number of nanostars in the systems is (d) N = 50, (e) N = 150, and (f) N = 250. (g–i) MSD of 5WJ nanostar mixtures. The number of nanostars in the systems is (g) N = 50, (h) N = 150, and (i) N = 250. (j) Diffusion constant of different nanostar systems at different temperatures and densities. The points are calculated using velocity autocorrelation and the line corresponds to the diffusion constant from the slope of the MSD vs. time graph.

IV. Conclusions

In summary, we have presented a bead-spring model of DNA nanostars to study the structure and governing dynamics and thermodynamics of their gelation. This simple model captures the essential qualitative feature of DNA gelation for a wide range of densities. Using RDF calculation, we have shown that all the systems transform from an unstructured fluid to a gel-like phase upon temperature lowering. The amorphous nature of the gel phase was confirmed using structure factor calculation. The phase transition from an arrested gel phase to an unstructured liquid phase has been modelled using a two-state theoretical model. We find that the loss of configuration and translational entropy is mainly giving rise to a gel phase at low temperatures. There is a competition between the entropic loss and the enthalpic gain in the gel phase. The absolute translational and rotational entropy of the system which is evaluated using a two-phase thermodynamic model, also substantiates the above speculation. The slowing down of the dynamics upon approaching the transition temperature demonstrates the phase transition to a gel phase. The present work offers a thorough understanding of DNA hydrogel formation for different complex nanostars. Also, in this work, we have made an effort to evaluate the absolute entropy of a model gel system and tried to explain the possible origin of gel formation from a thermodynamic perspective. We would like to mention that, in the present study, we tried to decipher the gelation process using the simplest model which can provide a reasonable description of the nanostars, and also at the same time capture the essential features of the gelation phenomena. The model is planar by definition. Thus, we do not introduce any angular flexibility into the model. There already exist many other coarse-grained models like oxDNA and Martini which have the flexibility and can be studied to see the effect of flexibility on gelation.

With regard to DNA hydrogels, our interest lies in the qualitative macroscopic structure, dynamics, and thermodynamics of large-scale systems while correctly describing the underlying microscopic mechanism. Here, we consider a class of materials, DNA hydrogels, that are self-assembled from multivalent building blocks of DNA nanostars. The coarse-grain model chosen is computationally efficient, straightforward, and robust enough to capture the self-assembly of gelation. Some simplifications have been made, like neglecting sequence specificity and electrostatic interaction. However, these can be incorporated by re-parameterization of the interaction strength. Even though there are more accurate CG models (oxDNA or Martini) that can be studied if required, which can meet these deficits, but at the cost of computational power. Overall, this model can be implemented to study the gel formation of similar other chemical species where the bond formation is mediated by hydrogen bonding. We believe that the study will be helpful in understanding the underlying physics of gel formation of hyper-branched DNA-like nanostars for various bio- and nano-technological applications.

Conflicts of interest

The authors have no conflict of interest to declare.

Acknowledgements

We thank TUE-CMS, IISc, Bangalore, funded by DST, for providing the CPU hours and DAE, India, for financial support. SN thanks IISc for RA fellowship. DB thanks SERB-DST for Ramanujan fellowship; DBT-EMR, Gujcost-DST, GSBTM for research grant and IITGN for initial startup support. SN thanks Jayeeta Chattopadhyay for critically reading the manuscript.

References

1. A. Aggarwal, S. Naskar, A. K. Sahoo, S. Mogurampelly, A. Garai and P. K. Maiti, Curr. Opin. Struct. Biol., 2020, 64, 42–50.
2. N. C. Seeman, J. Theor. Biol., 1982, 99, 237–247.
3. N. R. Kallenbach, R.-I. Ma and N. C. Seeman, Nature, 1983, 305, 829–831.
4. N. C. Seeman, Nature, 2003, 421, 427–431.
5. E. Winfree, F. Liu, L. A. Wenzler and N. C. Seeman, Nature, 1998, 394, 539–544.
6. P. W. Rothemund, Nature, 2006, 440, 297–302.
7. D. Bhatia, S. Arumugam, M. Nasilowski, H. Joshi, C. Wunder, V. Chambon, V. Prakash, C. Grazon, B. Nadal, P. K. Maiti, L. Johannes, B. Dubertret and Y. Krishnan, Nat. Nanotechnol., 2016, 11, 1112.
8. M. Langecker, V. Arnaut, T. G. Martin, J. List, S. Renner, M. Mayer, H. Dietz and F. C. Simmel, Science, 2012, 338, 932–936.
9. H. Lee, A. K. Lytton-Jean, Y. Chen, K. T. Love, A. I. Park, E. D. Karagiannis, A. Sehgal, W. Querbes, C. S. Zurenko and M. Jayaraman, Nat. Nanotechnol., 2012, 7, 389–393.
10. A. V. Pinheiro, D. Han, W. M. Shih and H. Yan, Nat. Nanotechnol., 2011, 6, 763.
11. N. C. Seeman, Mol. Biotechnol., 2007, 37, 246.
12. D. Bhatia, S. Surana, S. Chakraborty, S. P. Koushika and Y. Krishnan, Nat. Commun., 2011, 2, 1–8.
13. X. Dai, Q. Li, A. Aldalbahi, L. Wang, C. Fan and X. Liu, Nano Lett., 2020, 20, 5604–5615.
14. S. H. Um, J. B. Lee, N. Park, S. Y. Kwon, C. C. Umbach and D. Luo, Nat. Mater., 2006, 5, 797–801.
15. J. B. Lee, S. Peng, D. Yang, Y. H. Roh, H. Funabashi, N. Park, E. J. Rice, L. Chen, R. Long and M. Wu, et al. , Nat. Nanotechnol., 2012, 7, 816–820.
16. E. Cheng, Y. Xing, P. Chen, Y. Yang, Y. Sun, D. Zhou, L. Xu, Q. Fan and D. Liu, Angew. Chem., Int. Ed., 2009, 48, 7660–7663.
17. V. Morya, S. Walia, B. B. Mandal, C. Ghoroi and D. Bhatia, ACS Biomater. Sci. Eng., 2020, 6, 6021–6035.
18. Y. Shao, H. Jia, T. Cao and D. Liu, Acc. Chem. Res., 2017, 50, 659–668.
19. S. Walia, V. Morya, A. Gangrade, S. Naskar, A. Guduru Teja, S. Dalvi, P. K. Maiti, C. Ghoroi and D. Bhatia, ACS Biomater. Sci. Eng., 2021, 7, 5933–5942.
20. S. Nagahara and T. Matsuda, Polym. Gels Networks, 1996, 4, 111–127.
21. E. Cheng, Y. Xing, P. Chen, Y. Yang, Y. Sun, D. Zhou, L. Xu, Q. Fan and D. Liu, Angew. Chem., Int. Ed., 2009, 48, 7660–7663.
22. Y. Xing, E. Cheng, Y. Yang, P. Chen, T. Zhang, Y. Sun, Z. Yang and D. Liu, Adv. Mater., 2011, 23, 1117–1121.
23. J. B. Lee, S. Peng, D. Yang, Y. H. Roh, H. Funabashi, N. Park, E. J. Rice, L. Chen, R. Long and M. Wu, et al. , Nat. Nanotechnol., 2012, 7, 816–820.
24. C. Li, X. Zhou, Y. Shao, P. Chen, Y. Xing, Z. Yang, Z. Li and D. Liu, Mater. Chem. Front., 2017, 1, 654–659.
25. O. Okay, J. Polym. Sci., Part B: Polym. Phys., 2011, 49, 551–556.
26. C. Li, A. Faulkner-Jones, A. R. Dun, J. Jin, P. Chen, Y. Xing, Z. Yang, Z. Li, W. Shu and D. Liu, et al. , Angew. Chem., 2015, 127, 4029–4033.
27. C. Li, M. J. Rowland, Y. Shao, T. Cao, C. Chen, H. Jia, X. Zhou, Z. Yang, O. A. Scherman and D. Liu, Adv. Mater., 2015, 27, 3298–3304.
28. C. Li, P. Chen, Y. Shao, X. Zhou, Y. Wu, Z. Yang, Z. Li, T. Weil and D. Liu, Small, 2015, 11, 1138–1143.
29. J. Jin, Y. Xing, Y. Xi, X. Liu, T. Zhou, X. Ma, Z. Yang, S. Wang and D. Liu, Adv. Mater., 2013, 25, 4714–4717.
30. L. Rovigatti, F. Smallenburg, F. Romano and F. Sciortino, ACS Nano, 2014, 8, 3567–3574.
31. Y. Shao, H. Jia, T. Cao and D. Liu, Acc. Chem. Res., 2017, 50, 659–668.
32. X. Zhou, C. Li, Y. Shao, C. Chen, Z. Yang and D. Liu, Chem. Commun., 2016, 52, 10668–10671.
33. Z. Xing, A. Caciagli, T. Cao, I. Stoev, M. Zupkauskas, T. O'Neill, T. Wenzel, R. Lamboll, D. Liu and E. Eiser, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 8137–8142.
34. F. Bomboi, F. Romano, M. Leo, J. Fernandez-Castanon, R. Cerbino, T. Bellini, F. Bordi, P. Filetici and F. Sciortino, Nat. Commun., 2016, 7, 1–6.
35. S. Roldán-Vargas, F. Smallenburg, W. Kob and F. Sciortino, J. Chem. Phys., 2013, 139, 244910.
36. F. Romano and F. Sciortino, Phys. Rev. Lett., 2015, 114, 078104.
37. S. Biffi, R. Cerbino, G. Nava, F. Bomboi, F. Sciortino and T. Bellini, Soft Matter, 2015, 11, 3132–3138.
38. F. Bomboi, S. Biffi, R. Cerbino, T. Bellini, F. Bordi and F. Sciortino, Eur. Phys. J. E: Soft Matter Biol. Phys., 2015, 38, 1–8.
39. F. Bomboi, D. Caprara, J. Fernandez-Castanon and F. Sciortino, Nanoscale, 2019, 11, 9691–9697.
40. J. Fernandez-Castanon, F. Bomboi, L. Rovigatti, M. Zanatta, A. Paciaroni, L. Comez, L. Porcar, C. J. Jafta, G. C. Fadda and T. Bellini, et al. , J. Chem. Phys., 2016, 145, 084910.
41. J. Fernandez-Castanon, F. Bomboi and F. Sciortino, J. Chem. Phys., 2018, 148, 025103.
42. J. Fernandez-Castanon, S. Bianchi, F. Saglimbeni, R. Di Leonardo and F. Sciortino, Soft Matter, 2018, 14, 6431–6438.
43. J. Largo, F. W. Starr and F. Sciortino, Langmuir, 2007, 23, 5896–5905.
44. I. D. Stoev, T. Cao, A. Caciagli, J. Yu, C. Ness, R. Liu, R. Ghosh, T. ONeill, D. Liu and E. Eiser, Soft Matter, 2020, 16, 990–1001.
45. H. Joshi, A. Kaushik, N. C. Seeman and P. K. Maiti, ACS Nano, 2016, 10, 7780–7791.
46. S. Naskar, H. Joshi, B. Chakraborty, N. C. Seeman and P. K. Maiti, Nanoscale, 2019, 11, 14863–14878.
47. S. Naskar, M. Gosika, H. Joshi and P. K. Maiti, J. Phys. Chem. C, 2019, 123, 9461–9470.
48. S. Naskar, S. Saurabh, Y. H. Jang, Y. Lansac and P. K. Maiti, Soft Matter, 2020, 16, 634–641.
49. S. Naskar and P. K. Maiti, J. Mater. Chem. B, 2021, 9, 5102–5113.
50. T. E. Ouldridge, A. A. Louis and J. P. Doye, J. Chem. Phys., 2011, 134, 02B627.
51. P. Šulc, F. Romano, T. E. Ouldridge, L. Rovigatti, J. P. Doye and A. A. Louis, J. Chem. Phys., 2012, 137, 135101.
52. J. J. Uusitalo, H. I. Ingólfsson, P. Akhshi, D. P. Tieleman and S. J. Marrink, J. Chem. Theory Comput., 2015, 11, 3932–3945.
53. A. Savelyev and G. A. Papoian, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 20340–20345.
54. T. A. Knotts IV, N. Rathore, D. C. Schwartz and J. J. De Pablo, J. Chem. Phys., 2007, 126, 02B611.
55. C. Maffeo, T. T. Ngo, T. Ha and A. Aksimentiev, J. Chem. Theory Comput., 2014, 10, 2891–2896.
56. Z. Xing, C. Ness, D. Frenkel and E. Eiser, Macromolecules, 2019, 52, 504–512.
57. W. Humphrey, A. Dalke and K. Schulten, et al. , J. Mol. Graphics, 1996, 14, 33–38.
58. B. B. Laird and A. Haymet, J. Chem. Phys., 1992, 97, 2153–2155.
59. S.-T. Lin, M. Blanco and W. A. Goddard III, J. Chem. Phys., 2003, 119, 11792–11805.
60. S.-T. Lin, P. K. Maiti and W. A. Goddard III, J. Phys. Chem. B, 2010, 114, 8191–8198.
61. T. A. Pascal, S.-T. Lin and W. A. Goddard III, Phys. Chem. Chem. Phys., 2011, 13, 169–181.

---

Footnote

† Present address: Institut Charles Gerhardt Montpellier (ICGM), University of Montpellier, CNRS, ENSCM, 34095 Montpellier, France.

---