Among the millions of biochemical processes and physicochemical activities involved in the proper functionality of living matter, perhaps one of the most interesting phenomena is the structural organization of double-stranded DNA (dsDNA) chains into some specific type of bio-macromolecules [1–6]. All of them protecting in some manner the genome from external agents. From the atomic composition of dsDNA, it is well-known that this polyelectrolyte has an effective negative charge due to the phosphate groups along the chain [7]. This type of charge distribution causes that the DNA interacts electrostatically with itself and other molecules, such as proteins, enzymes and ions, leading to a complex self-organization of the DNA chains that is linked to their functionality. Thus, the full understanding of the DNA properties under controlled experimental conditions is critical to improve our knowledge of its function and relation to certain biochemical activities.

Under bulk conditions, the pioneering contributions of Rau and Parsegian, who carried out X-ray diffraction experiments for long (length ≫ cross section) and parallel DNA strands, elucidated the effects of osmotic pressure on the interstrand distance and the emergence of the so-called hydration forces [8, 9], as well as the consequence of varying the ionic conditions on the interstrand interactions, which can be switched from repulsive to attractive to promote, for example, DNA condensation, depending on the polyvalency of the buffer [10, 11]. These experiments were carried out under thermodynamic conditions that allowed the interchange of water and ion molecules with the buffer in order to keep an homogeneous chemical potential at several polyethylene glycol concentrations, which established an osmotic pressure that induced the ordering of DNA strands.

In Biology, there are some examples of bio-particles, such as cells, sperm and viruses that, in addition to their complex internal components, have DNA chains with a characteristic length [12]. Notably, this DNA is tightly packed into a small volume, which gives rise to a high DNA packing density, whose structural arrangement is influenced by the temperature, interstrand interaction (mediated by the ionic strength of the buffer) and DNA persistence length (bending stress). Even more, the internal pressure of such biomacromolecules plays a key role in the appropriate function of them as the case, for example, of the lytic life cycle of viruses [13]. For example, in the wild-type bacteriophage lambda and herpes simplex virus type 1 [14, 15], due to the high packing of the viral DNA into the capsid, there is experimental evidence of the huge internal pressure, which is of some tens of atmospheres. This high pressure environment facilitates the rapid ejection of the DNA from the capsid into the host cell, which will be infected to initiate the process of viral replication.

Providing knowledge on the physical properties of DNA at high pressures, either hydrostatic or osmotic, could help in, for example, the development of new technologies, such as antiviral molecules [16], as well as to explain pathological processes that might be associated with cancer [17, 18] and the formation of new vectors for gene therapy [19, 20]. Another highlighting phenomenon is the so-called pressure-driven dsDNA denaturation process, which is mainly related to ionic strength variations [21, 22]. In this case, the increase of ionic concentration causes an osmotic pressure difference, which corresponds to a completely different physical scenario than the Rau and Parsegian’s experiments discussed above. In particular, pressure-driven dsDNA denaturation has been described in terms of a thermal-electrostatic competition, i.e., under certain thermodynamic conditions, thermo-mechanical variations can lead to DNA denaturation [21, 22]; pressure-induced DNA melting is basically achieved when the breaking of the hydrogen bonds between base pairs takes place. This mechanism can also be associated with irreversible damage caused to DNA that might be linked to the onset of cancer [23]. This issue plays a key role in the transcription of genomic information, particularly in the synthesis of other biomolecules.

On the other hand, in recent years, there has been a growing interest in the self-assembly and the thermodynamic behavior of bioparticles, such as lysozyme [24–26] and monoclonal antibodies [27, 28], to investigate the mechanisms that give rise to, for instance, the phase separation and aggregation (or particle clustering); processes that are crucial and may affect the functionality of such biomolecules [29, 30]. Typically, the interaction between bioparticles is mediated by the heterogeneous surface charge distribution [31–33], i.e., proteins and enzymes have regions or specific active sites that are electrically charged to recognize other agents in order to bind them, and lock them or catalyze them, mainly, to complete a specific biochemical function [34–36]. This type of process implies a specific orientation or spatial organization of the bioparticles in the available volume, which should also depend on the thermodynamic conditions. This contribution explores the organization and stability of two double-stranded and electrostatically interacting DNA chains under high hydrostatic pressures, i.e., our interest is focused on the hydrostatic pressure exerted by solvent molecules over DNA chains and its effect on their spatial organization. To this end, we have performed all-atomistic molecular dynamics (MD) simulations of two short DNA chains to provide some preliminary insights into the effects of high pressurization and finite particle size on the self-assembly of DNA chains. The MD simulations include explicitly the presence of water molecules and ions to emulate, as close as possible, the physiological conditions where the DNA is typically embedded. We have also explored hydrostatic pressures that are biologically relevant.

The results discussed in this work are based on atomistic MD trajectories for two interacting double-stranded B-DNA dodecamers using the parmbsc1 force field [37–39] in conjunction with the GROMACS software package version 2019 [40]; parmbsc1 can be found in GROMACS as amber99bsc1. The parmbsc1 is a general-purpose atomistic force field that takes into account high-level quantum mechanical data [37] and has been tested over a wide variety of DNA sequences. In a coarse-grained picture, we have simulated two strongly screened like-charged flexible rods, immersed in a salty aqueous medium, whose interaction at the cytosine and guanine ends is of electrostatic nature.

The sequence and starting structure of each oligomer were obtained from the Protein Data Bank with code 1BNA. This DNA structure is typically used as a benchmark for the development of force-fields [37]. Each MD run was performed using two double-stranded B-DNA dodecamers (d(CGCGAATTCGCG)2) initially separated by 3 nm and centered in a cubic box of around 8.5 nm, i.e., 2.1 the linear size of a DNA chain (see Figure 1). Although the cubic simulation box has the computational disadvantage of requiring a greater number of water molecules, as compared to other simulation cells with different geometrical shape, it gives rise to the same spatial freedom in the three spatial directions. Therefore, the dodecamers have the same possibility of translational and rotational motion. In addition, other simulation boxes, such as the rhombic dodecahedron and truncated octahedron, are more suitable for globular biomolecules. Several initial conditions and starting DNA arrangements, i.e., parallel, perpendicular and oblique, were used in order to disregard any effect associated with the initial configuration. Structures were then solvated with the TIP3P water model with a minimum of 1.2 nm buffer solvation layer beyond the solute, neutralized, and 0.5 M excess NaCl was added. A typical simulated system consisted of around 18600 water molecules and a total number of 58000 atoms. Larger systems were also considered to discard finite size effects (data not shown). DNA chains, water and salt molecules are in a closed thermodynamic system, so we are essentially testing the effects of hydrostatic pressure.

MD simulations using the standard leap-frog integrator at the absolute temperature of 300 K and six different hydrostatic pressure values, namely, p = 1, 10, 20, 30, 40, and 50 bar, were carried out using the Verlet cut-off scheme [40]; periodic boundary conditions were applied in all spatial directions. An integration step of 2 fs in conjunction with SHAKE to constrain and fix the X-H bonds to the corresponding default values were considered [41]. The Lennard-Jones and Coulomb interactions were used and computed within a cutoff radius of 1.0 nm, and long-range electrostatic interactions were calculated using the particle mesh Ewald (PME) method [42] with cubic interpolation; GROMACS automatically tunes the load balance by scaling the short-range electrostatic cutoff and grid spacing.

The water-ion-DNA system at 1 bar was optimized using the standard steepest descent energy minimization, followed by a thermalization process during 5 ns using the NVT canonical ensemble. During this equilibration, the temperature was kept constant using the velocity-rescale thermostat [43], with the restriction of fixing the position of heavy atoms to equilibrate the solvent and ions around both DNA molecules. The last configuration of this process was used as the initial configuration in the following steps. After reaching the desired equilibrium temperature, a further pre-equilibrated run for 5 ns in the NPT isobaric-isothermal ensemble was performed using the so-called Berendsen barostat [44] to reach in each case the equilibrium hydrostatic pressure. Subsequently, we have equilibrated each water-ion-DNA system at the targeted hydrostatic pressure for an additional 10 ns period using the Parrinello-Rahman [45] barostat with a 2 ps relaxation time and compressibility of 4.5 × 10⁻⁵ bar⁻¹. Finally, a run of 240 ns for each water-ion-DNA system, i.e., p = 1–50 bar, was carried out for gathering statistics and production. We make sure that both DNA chains freely rotate their own size in a timescale shorter than the time window used in the MD simulations. Conformational configurations were saved every 50 ps given a total of 4800 configurations, and visualization snapshots were prepared with the VMD software [46].

We now show the results of the MD computer simulations of two short DNA chains immersed in an explicit aqueous solution under different hydrostatic pressures. This work focuses on the evolution of the relative separation between DNA chains. This information is explicitly reported in Figure 2. We show the results that were obtained by considering that both DNA molecules are initially distributed in a parallel configuration and located at a relative distance (measured from the centers of mass (COM)) of 3 nm. However, as mentioned above, we have used other initial configurations and considered larger simulation boxes (keeping fixed the total density) within an extended time window, obtaining practically the same result (data not shown). As seen in the figure, in most pressure conditions, both chains reach an equilibrium co-linear configuration, except for the p = 10 bar case (see snapshots at t = 240 ns). The fact that the two dsDNA chains align along the same main axis is, of course, a consequence of the finite size of the chains. This phenomenology contrasts with the completely different scenario reported for the DNA at the bulk under high osmotic pressure conditions [8, 9], where the osmolyte induced a parallel ordering of the DNA strands. Although hydrostatic and osmotic pressures are not the same, either type can induce a chain-like arrangement of DNA oligomers. Thus, the structural evolution towards the co-linear state clearly depends on the thermodynamic conditions, i.e., hydrostatic pressure. No mechanical instability or denaturation of the DNA molecules is observed. Interestingly, at low and high pressures, the relative distance between chains is the same, about 4 nm, whereas at p = 10 bar, it reaches a value of 3 nm.

More explicitly, at p = 1 bar (see Figure 2A), the structural reorganization of the DNA molecules occurs monotonically and very fast (at around 10 ns), and the fluctuations around the equilibrium state do not exceed 0.25 nm. This behavior can be understood as follows. Under this pressure condition, the couple of dsDNA molecules start to evolve from a quasi-parallel configuration; however, the electrostatic interactions strongly dominate over the solvent-mediated interaction, which causes the fast rearrangement of chains, looking for the most favorable configuration that minimizes the energy; see Figure 3A. On the other hand, at p = 10 bar (see Figure 2B), the evolution towards the co-linear configuration is not monotonous anymore and, in fact, it is not the predominant one, but, instead, an L-type conformation appears. For this particular hydrostatic pressure condition, one can clearly appreciate that the DNA biomolecules explore an extended configurational landscape and it seems that during an extended time window, they try to be far away from each other. This scenario can be linked to the competition between the electrostatics and the solvent-mediated interaction, i.e., mechanical stresses; the latter enhanced by the pressure; see Figure 3B. Furthermore, the structural relaxation takes a much longer time ( ∼ 80  ns) as compared with the previous case, but now the DNA chains are closer each other although the relative distance between them exhibits stronger variations around the equilibrium state; fluctuations that occur at a distance of roughly 3 nm measured from the COM. As seen in Figure 3, at all hydrostatic pressure conditions, the electrostatic contribution clearly dominates over the other ones.

It is quite interesting the fact that an increase of one order of magnitude of the pressure value induces a slower structural relaxation. Nevertheless, if the pressure is further increased, Figures 2C, D, which corresponds to pressure values of p = 30 and 50 bar, respectively, the chains are twisted to each other and, at earlier times, they try to explore, again, a rich set of configurations in order to reach the most stable one. The latter being, once again, a co-linear spatial organization whose chain distance fluctuations are much smaller than the previous pressure conditions, describing a persistent equilibrium conformation; see, e.g., Figure 3C. We should point out that the results presented in Figure 2 are very robust in the sense that different initial conditions and several starting DNA configurations were tested and the outcome in all cases was basically the same, even the case p = 10 bar exhibited the same L-type configuration regardless the initial configuration or distribution of dsDNA molecules. In appendix, we discuss in more detail the effect of the initial conditions and the biomolecular recognition on the final configuration state reached by the dsDNA molecules at p = 10 bar.

The values of the hydrostatic pressure covered in the MD simulations and explored in this work are commonly found in biological environments; larger pressure values, although of academic interest, might not be realistic. Under these conditions, since the density of the whole water-ion-DNA system is basically determined by the water density, as expected, the volume (or dimensions of the simulation box) where the dsDNA molecules are immersed did not change significantly (data not shown), i.e., the water behaves as an incompressible fluid under the thermodynamic conditions here considered. However, it was interesting to note that the relaxation of the internal stresses (produced by the hydrostatic pressure) in the fluid clearly affected the exploration of the configurations that minimize the energy. As discussed in the previous section, except for p = 10 bar, the equilibrium state in the other cases corresponded to a co-linear configuration, which resulted from a delicate interplay of electrostatic interactions and hydrogen bonding between base pairs; the kinetics or speed of the exploration of the configurational states seemed to be mainly determined by both the hydrostatic pressure and molecular recognition, as discussed further below. Other cases at p = 20 and 40 bar were also studied and the co-linear configuration was again the equilibrium state achieved by the dsDNA chains (data not shown). This reflects the generality of the preliminary results reported in this contribution. We should also point out that the fluctuations of the COM could be further elucidated by characterizing the relative dynamics of the two DNA dodecamers, i.e., it is appropriate to determine if those fluctuations are induced, for instance, by hydrodynamic correlations or if they are entirely inherent to the DNA backbone of each dodecamer. This point will be reported elsewhere.

More importantly, we have observed a co-linear arrangement during the different pressure conditions with a non-specific interaction. For example, in the situations at p = 1 and 30 bar, the rearrangement between the DNA dodecamers is generated between terminal cytosine-guanine chains A and B, respectively, but at p = 50 bar it shows that the base pair interaction was between terminal cytosine-cytosine and guanine-guanine of the A and B chains, respectively. The most interesting case was at p = 10 bar, where the pair base interaction is, as already discussed, like an L form. This kind of interaction is due to the structural complementarity or biomolecular recognition, which seemed to result from the pressure exerted on the water-ion-DNA system. However, as displayed in Figure 3, the electrostatic energy became the most important energetic contribution of the water-ion-DNA system, then another aspect that could shed light on the feasibility of the L-shape configuration observed for p = 10 bar is the spatial organization of the ions around the dodecamers. In Figure 4, two snapshots of the water-ion-DNA system are displayed. As illustrated in the figure, the Na⁺ ions are preferentially distributed on the surface of the DNA chains, indicating that they dominate on the DNA-DNA interaction, as expected. The L-type conformation in dsDNA chains has also been observed in DNA molecules confined between lipid membranes [47]. We should stress out that, however, at all pressure conditions here studied, intrastrand complementarity (G-C) was preserved at the sites of co-linearity of the dsDNAs (see snapshots of Figure 2).

During the analysis of the structural evolution of both dsDNA molecules, we also focused on the local structural changes of the typical conformation of B- to non-B-DNA, since it is well-known that they can be induced by changes in the environmental conditions, protein effects, and superhelical tension [48]. These conformational variations directly affect replication, gene expression, recombination, and mutagenesis. Our results indicate that the dodecamers maintained a DNA canonical structure, keeping backbone angles and helical parameters within the expected values for a B-DNA duplex (data not shown). In general, highly charged polyelectrolytes, such as duplex DNA, exhibit significant chain stiffness compared to neutral polymers of similar molecular length due to electrostatic repulsion arising from localized point charges along the molecular boundary. Within our preliminary findings, we have observed a different behavior at p = 10 bar. Under this pressure condition, the water-ion-DNA system reached an equilibrium configuration at a similar time to the other states, but the relative distance between the COM of the dodecamers is 1 nm lower, as a result of the arrangement in L between the dsDNA molecules.

On the other hand, the sugar-phosphate backbone of DNA is hydrophilic. The internal portion of the bases (the rings) has less polar characteristics, although not totally hydrophobic, than the rest of the molecule [49]. The interactions generated between the faces of the rings of the bases of each nucleotide are of the “stacking” type and are highly stable [50, 51]. These aromatic interactions have been proposed to consist of vdW scattering, electrostatic, and hydrophobic forces, but each component’s relative contribution and magnitude still need to be determined [50, 52]. There are three reported representative conformations of pi-stacking interactions with different energy of conformation. These conformations include face-to-face (−1.48 kcal/mol), edge-to-face or T-shaped (−2.46 kcal/mol) and offset (−2.48 kcal/mol) [52]. The offset interaction is present in the interior of the dsDNA chain and on the edge of the co-linear interactions observed in our MD simulations, since it represents the lowest conformation energy. In contrast, a ring-interacting T-shaped conformation seems to be responsible of the L-shaped conformation at p = 10 bar. In this sense, the co-linear conformation generated by the interaction between two ends in the dsDNA sequences is more stable than that created by a single interaction, such as in the L-shaped conformation. Therefore, our hypothesis is that the L-shaped conformation corresponds with a local minimum that is very close to the global one (see appendix). This hypothesis will be tested with further computer simulations.

Although the outcome of this contribution is centered on the very dilute regime at biologically relevant hydrostatic pressures, in contrast, there have been others all-atomistic computer simulation efforts to understand the phase transition of DNA molecules at finite concentration. Most notably, a novel computational method on a multiscale simulation based on AdResS [53, 54] has shown some interesting results at different DNA densities. In particular, it has been reported the influence of the osmotic pressure on the transition of the DNA spatial organization from cholesteric to hexagonal to orthorhombic arrays [55–57], as well as the effect of the ionic strength, concluding that hydration forces are the dominant interaction at high DNA packing, as experimentally highlighted a while ago [10, 11]. This confirms the importance of the explicit inclusion of solvent molecules in order to consider, as much as possible, all fundamental interactions involved in the phase state of any kind of biomacromolecule. Even though the high stress on DNA chains via hydrostatic or osmotic pressures describes different biophysical scenarios, both of them clearly influence the DNA arrangement into the available volume and shed light on our understanding of the thermodynamic conditions that the DNA faces into living matter.

As stated above, we have reported preliminary results that point toward an interesting structural evolution of interacting dsDNA molecules under the influence of hydrostatic pressure. Of course, there are many routes that can be further explored. For example, to carry out a systematic analysis of the pressure-induced self-assembly of a larger number of dsDNA chains, to study longer DNA sequences to mimic the biological conditions and to examine the effect of hydration forces on the structural organization of dsDNA. Work along these lines is currently in progress.