Molecular dynamics simulation of complexation between plasmid DNA and cationic peptides

The elucidation of the process by which cationic peptides condense plasmid DNA (pDNA) is important for unraveling the mechanism of peptide/pDNA complex formation, which plays a vital role in gene delivery for the genetic transformation of living cells. We performed atomic MD simulations of the complexation of pDNA in the presence of two cationic peptides, KH9 (with an alternating sequence of lysine and histidine) and Cytcox (functioning as a mitochondria-targeting signal), to investigate the mechanism of pDNA condensation. The simulations revealed that the cationic peptides bound to the pDNA and that defects in pDNA formed in response to the densely packed cationic peptides, presumably initiating the folding of the double-stranded pDNA into a globule morphology. The decrease in the radius of gyration and the number of hydrogen bonds and the increase in the writhe structure, with a slightly higher tendency for the Cytcox/pDNA system, strongly support the formation of pDNA defects leading to the bending of the double helix. The results provided insight into the mechanism of pDNA complexation with cationic peptides, which should contribute to the future design of highly efficient gene delivery systems using peptide-mediated nanocarriers. Atomic molecular dynamics simulations of the complexation of plasmid DNA in the presence of two types of cationic peptides clear the mechanism of plasmid DNA condensation for gene delivery systems.


Introduction
Deoxyribonucleic acid (DNA) is a biopolymer responsible for carrying genetic information and is organized into a rigid double helix by interactions between multiple base pairs of two polynucleotide chains.The anionic charge of each nucleotide unit allows DNA to form complexes with cationic counterparts, such as histones, which pack DNA into chromatin in nature [1].The complexation of DNA with cationic polymers, including cationic peptides, has been extensively studied for gene delivery applications in living cells [2].The electrostatic interaction between DNA and cationic peptides suppresses the electrostatic repulsion of anionic nucleotides, leading to the formation of complexes suitable for gene delivery with size compaction.We have previously developed cationic peptide-based nanocarriers, which consist of a cationic peptide and a functional peptide such as a cell-penetrating peptide or organelle targeting signal, to deliver plasmid DNA (pDNA) into plant cells [3][4][5].The DNA/peptide complex formed a spherical assembly with a hydrodynamic size ranging from tens to hundreds of nanometers.The size of the complex is a critical factor influencing the internalization process through cell walls/cell membranes to increase the efficiency of gene delivery into the cells [6].Therefore, unveiling how cationic peptides condense DNA is of interest to facilitate the design of peptides that produce homogeneous condensed particles for gene delivery applications and could contribute to further development for efficient genetic modification by gene delivery systems.
The double-stranded DNA is rigid with a long persistent length of approximately 50 nm [7].Interaction with cationic peptides partially unwinds the double helix into single strands, which causes the transition from coil to globule.Bending of the double helix should occur during the condensation process, although analytical methods that can directly investigate the DNA condensation process by cationic polymers are quite limited.The static morphology achieved by pDNA condensation by cationic peptides has been visualized by microscopic observation using fluorescence microscopy or atomic force microscopy (AFM) [7][8][9].However, DNA circles of plasmid size have been proven to be too structurally disordered to be crystalized for X-ray crystallography and too large to visualize by NMR.The condensation process of the pDNA/peptide complex has yet to be fully clarified.In this context, molecular dynamics (MD) simulation of DNA molecules is a suitable technique to follow the morphological changes during the condensation process at a molecular level [10].Small DNA minicircles, due to their relatively small size, have proven to be a good model to study DNA supercoiling using a combination of experimental techniques and MD simulations using implicit and explicit solvents [10].In this study, we performed MD simulations of the complexation of pDNA with cationic peptides to visualize the process of DNA condensation.

Materials and methods
The creation of circular DNA and the protocol for molecular dynamics simulation in implicit solvent were performed following the protocol of Noy et al. [10].The Ambertools program Nucleic Acid Builder (NAB) was used to create an atomistic model of a circular DNA with 260 base pairs (bp) and Δlinkingnumber (supercoiling) of -1, which is defined as the number of extra turns to add or subtract to the ideal B-DNA conformation of 10 base per turn.The obtained pdb file was used to construct the DNA system with Ambertools tleap [11].
For the simulation of peptides and DNA, the combination of the amber ff14SB and bsc1 force field was used since it has been reported to give the most useful results [10].A bond between the 5ʹ and 3ʹ termini of the DNA was manually created using tleap to obtain a closed circular DNA structure, and we created the topology and coordinate files needed to start MD simulations with AMBER.
Energy minimization and equilibration were performed in Amber18 using the protocol of Noy et al. [10].The default PBradii was set to 'mbondi3', which is the recommended radii for the intended generalized Born (GB) solvation model to use for protein and DNA simulation [10] with the Langevin thermostat for the stages of the simulations with a nonbonded cutoff set to 60 Å at 300 K and 200 mM salt concentration [12].The parameter "gam-ma_ln" was set to 0.01 ps −1 to reduce the collision frequency calculated for water and the reduced water viscosity [10].The systems were equilibrated in a series of stages.
First, energy minimization was performed for 10000 cycles with all atom restraints of 100 kcal mol −1 Å −2, followed by an additional 10000 cycles without restraints.After that, a three-stage equilibration was performed on the systems with all-atom restraints.In the first stage, the system was linearly heated for 400 ps from 10 to 300 K with all atoms restrained by 50 kcal mol −1 Å −2 .Then, the system was run at a constant 300 K with reduced restraints for 100 ps, first with 10 kcal mol −1 Å −2 restraint and later with 1 kcal mol −1 Å −2 .
Systems of DNA with peptides were constructed from a simulated DNA system after 30 ns extracted from a clustered analysis.The procedure to create the systems with peptides was similar to the one described previously for DNA, with the difference that we used the protein force field ff14SB to describe the incorporated peptides.The NMR structure of Cytcox-KH9 (MLSLRPSIRFFK-KHKHKHKHKHKHKHKHKH) was used as the starting structure for creating both starting structures of KH9 and cytcox [13].The peptides were randomly located around the DNA at the calculated N/P ratios, defined as the number of charged amino acid residues (K and R) of the peptides to the number of phosphate groups of the pDNA.Thus, for the DNA-KH9 system at N/P = 0.5, 29 KH9 peptides were incorporated, and 87 peptides were incorporated into the DNA-Cytcox system at N/P = 0.5.Systems were energy minimized and equilibrated using the same protocol for the DNA.
The simulations were run for 100 ns after equilibration in implicit solvent using AMBER18.Implicitly solvated MD simulations were performed using GROMACS version 2020 with standard MD protocols for 100 ns [14].
Analysis of the simulations was carried out using CPPTRAJ analysis simulation software including Ambertools [15], and the characteristic parameters such as the radius of gyration (R g ) were determined.The geometrical descriptions of the writhe and the bend structures within the pDNA molecule were obtained using the WrLINE molecular contour and SerraLINE programs [16].The bending angles were calculated in the last 20 ns simulations using two tangent vectors separated by 16 nucleotides, a compromise length for capturing the overall bend produced by a defect [10].

Results
We previously designed and used cationic fusion peptides to assemble DNAs into spherical assemblies for gene delivery into plant, animal, and bacterial cells [2-5, 17].One component of the fusion peptide is an alternating sequence of lysine and histidine residues (KH9) that effectively condenses pDNA of interest into the spherical assembly.The other part of the fusion peptides are functional sequences, such as cell-penetrating peptides and organelle-targeting peptides, depending on the purpose of delivery.In this study, we focused on DNA condensation by each fragment of the fusion peptide to track the condensation process by MD simulation.As a functional fragment, we used a mitochondrial targeting sequence Cytcox, a partial presequence of yeast cytochrome c oxidase subunit IV, for DNA complexation [13].
The pDNA model with a relatively small size of 260 bp was constructed for efficient simulations at the atomistic level (Fig. 1a) [10] because the size of the pDNA previously used in our gene delivery studies was ~10 kbp, which is beyond the reach of atomistic MD simulations [3-5].The MD simulations in implicit solvent showed that even at short timescales, the DNA minicircles underwent fast structural rearrangements.This was because the conformational fluctuations of the systems were accelerated by at least one order of magnitude in generalized Born continuum solvent due to the absence of viscous damping, allowing a rapid sampling of different configurations of the pDNA [10].We analyzed the structure of the 3 systems, pDNA without peptides, KH9/pDNA, and Cytcox/pDNA, modeled with MD simulations in implicit solvent using a conformational clustering algorithm to select representative structures of the DNA [15].The KH9/pDNA and especially the Cytcox/pDNA had representative cluster conformations that appeared to be more compact and twisted than the pDNA only system (Fig. 1).The pDNA-only system mostly maintained the initial circular morphology throughout the simulation time, with the molecular plane slightly bent (Fig. 1d).When the local helical structure was examined, the double helix was slightly distorted, but almost all the base pairs maintained the hydrogen-bonded structure (Fig. 1g).On the other hand, in the pDNA system accommodating the KH9 peptides, the peptide molecules were strongly bound to the helical backbone throughout the pDNA by electrostatic interactions (Fig. 1e).Because the KH9 peptides were densely packed with nucleotide moieties, severe distortions and bending of the helical backbone occurred sporadically at some sites.The peptide molecules interfered with the proper hydrogen bonding of the base pairs at the locally distorted sites (Fig. 1h).The results indicated that the system reasonably simulated the initial stage of pDNA condensation by folding of the DNA backbone induced by interaction with the cationic peptides.
The pDNA system with the functional peptide Cytcox exhibited a similar trend to the KH peptide [5].Most of the peptides bound tightly to the pDNA backbone even with fewer cationic residues (Fig. 1f).The pDNA molecule was wholly packed with the peptides due to the larger number of Cytcox peptides.The distorted structures were also found sporadically throughout the pDNA molecule with hydrogen bonding of base pairs interrupted with peptides (Fig. 1i).The Cytcox peptide, which functions as a mitochondrial targeting signal in plant cells, contains three cationic residues, causing interaction with pDNA during complexation [13].In summary, the representative cluster conformations of KH9/DNA and especially Cytcox/DNA, appeared to be more compact and twisted than the DNA only system (Fig. 1).We also performed MD simulations of peptide/pDNA with different pDNA sizes of 100 or 340 bp and confirmed a similar complexation process.Our previous study revealed that the interaction with pDNA disrupted the original function of the functional segment (containing cationic residues) of the fused peptide nanocarriers, leading to the deterioration of gene expression efficiency in plant cells [18].Therefore, the simulation result with Cytcox indicated that appropriate design of the peptides is needed for multifunctional nanocarriers that condense DNA for gene delivery systems.
To quantitatively assess the DNA condensation ability of the cationic peptides, we analyzed various morphological factors regarding the global conformational change, including the number of writhes (coiling deformation of the helical backbone) [19], using the WrLINE program that extracts the molecular contour of DNA to calculate topological writhes [16].The smaller size of KH9/pDNA and Cytcox/pDNA was confirmed by the calculated radius of gyration of the DNA structures over the simulation time (Fig. 2).Compared with the pDNA-only system, complexation with the peptides tended to decrease the R g of the pDNA (123.3 ± 5.3 Å), with the Cytcox system having the smallest R g (102.6 ± 4.3 Å) (Figs. 2a, d).This indicates that the electrostatic interaction between the polynucleotide backbone and the peptides triggered the initiation of DNA condensation.The writhe structure is generally observed in supercoiled DNAs, which compensates for steric constraints of supercoiling.The pDNA system with Cytcox showed a greater propensity to form the writhe structure (0.79 turns) even with a small pDNA size than naked DNA (0.27 turns) and KH9/pDNA (0.39 turns).This is ascribed to the high constraint within the pDNA backbone caused by the large number of Cytcox molecules.While both peptides cause some degree of compaction on the DNA structure, at the same N/P ratio used here (0.5), Cytcox induces larger DNA compaction, as evidenced by the lower radius of gyration and writhe number, probably by the higher number of peptides (87) than KH9 (29).The number of intramolecular hydrogen bonds in the DNA molecule, associated with the hydrogen bonding of base pairs, decreased over the simulation time for the systems with KH9 and Cytcox (Figs. 2c and 2f).Both KH9 and Cytcox peptides disrupted the hydrogen bonding of base pairs, with a mean of 427.8 and 432.8 hydrogen bonds, respectively, compared to the naked DNA with a mean of 448.8 hydrogen bonds for all simulation times (Fig. 2c, f).The disruption of hydrogen bonding is thought to be caused by the interaction with cationic residues of the peptides.Therefore, the degree of hydrogen bond disruption almost comparable between KH9 and Cytcox due to the same N/P ratio of 0.5.results in the distortion of the double structure of pDNA observed in Fig. 1.The disruption to hydrogen bonding collapsed the complementary base pairing in the double helix, resulting in distorted, bent structures.
The bending angles of the pDNA backbone were calculated from the WrLINE profiles using the SerraLINE program.The pDNA-only system, where the circular morphology was maintained throughout the simulation, showed average bending angles ranging from 10°to 40°, whereas large bending angles up to 100°were sporadically obtained for the systems with KH9 and Cytcox (Fig. 3).It was reported that the bending angle is highly correlated with defect formation in pDNA with a critical threshold angle [10].The defects of pDNA in the MD simulation can be defined as a disruption to base stacking and complementary base pairing.Large bending angles of the DNA did not necessarily indicate the formation of a defect, but some of them matched the defects observed in the simulated conformations, and some defects appeared and then after some time returned to a normal conformation.
The morphological analyses of the MD simulations revealed that complexation with cationic peptides caused conformational changes disrupting the double helix structure.The condensation behavior of pDNA with cationic polypeptide-based block copolymers was experimentally examined in a previous study.Using shorter polypeptide segments, the pDNA underwent a regularly folded conformation with dissociation of the double helix only at each fold [20], and increasing the polypeptide length induced a morphological change to a globule structure [7].Our results using small pDNA indicate that the MD for pDNA with cationic peptides simulated the initial stage of folding of double-stranded pDNA that leads to a more condensed morphology.When we used fusion peptides to condense a large pDNA, typically approximately 10 kbp in size, the size of the peptide/pDNA complex was tens to several hundred nm [3-6, 13, 18].The MD simulation in this study also showed initial condensation behavior for peptide/ pDNA systems, although the size compaction is limited due to the rigidity attributed to the smaller pDNA size (260 bp).

Summary
We performed atomic MD simulations of the complexation of pDNA in the presence of cationic peptides to investigate the mechanism of pDNA condensation.The simulations revealed that the folding of the double-stranded pDNA into a globule morphology is presumably initiated by the formation of defects in response to the densely packed cationic peptides, leading to the bending of the double helix.Not only cationic peptides but also functional peptides containing cationic residues induced the bent structure.The results provided insight into the mechanism of pDNA complexation with cationic peptides, which should contribute to the future design of highly efficient gene delivery systems using peptide-mediated nanocarriers.

Fig. 1
Fig.1The 260 bp pDNA constructed in this study (a).Initial systems of the equilibrated pDNA with KH9 (b) and Cytcox (c) peptides.Snapshots of the simulations (100 ns) for the pDNA systems without peptides (d) and with KH9 (e) and Cytcox (f).The expanded images in the dashed squares are shown in (g-i)

Fig. 2
Fig. 2 Conformational analysis of DNA showing the calculated radius of gyrations (R g ) (a), number of writhes (b), and number of hydrogen bonds (c) for the pDNA systems with and without KH9 or Cytcox