Computer simulation of complex of lysine dendrigraft with molecules of therapeutic KED peptide

Lysine dendrimers and dendrigrafts are often used in biomedicine for drug and gene delivery to different target organs or cells. In present paper the possibility of complex formation by lysine dendrigraft and 16 molecules of therapeutic KED peptide was investigated using molecular dynamics simulation method. A system containing of one dendrigraftt and 16 KED peptides in water were studied. It was shown that stable complex consisting of the dendrigraft and the peptide molecules formed and structure of this complex was studied. Similar complexes could be used in future for delivery of other therapeutic peptides to


Introduction
Interest to macromolecules with regular brunched structuregrows every year [1]. One of the most wellknown polymers with dendritic structures are dendrimers. The use of different types of dendrimers in drug delivery research is discussed in many papers [2,3]. Dendrigraft is also a brunched polymer. Dendrigrafts could be described from one hand as dendrimers with short linear chain in their core or from another hand as dendritic brush with short main chain and long side chains. Lysine dendrigrafts consists entirely of lysine aminoacid residues (that are biocompatible) [4,5]. At the same time their terminal groups could be functionalized by other aminoacids or by other active groups or molecules [6]. Lysine dendrigrafts are polymers that are rich with amines. Due to this reason they could be used as antibacterial [7] or antiviral agents [8]. Also they could make complexes with oppositely charged peptides due to strong electrostatic interaction between their positively charged groups (NH3+) and negatively charged aminoacid side groups (COO-) of peptides. Hydrogen bonds between dendrigraft and peptide and hydrophobic interactions between their nonpolar groups are also important for creation of such complexes. Due to these ability to make complexes thedendrigrafts like other dendritic molecules could be used as.
Therapeutic KED peptide (Lys-Glu-Arg) [10][11] was selected for our study as a model peptide. This peptide belongs to a class of regulatory peptides and has an protective properties.
The goal of this paper is to demonstrate the possibility of complex formation between lysine dendrigraft of second generation and therapeutic KED peptide using molecular dynamics simulation method.

Molecular dynamics method
Molecular dynamics (MD) method is currently the main method for simulation of polymer and biopolymer systems. The method consists in numerical solution of the classical Newton equations of motion for all atoms of the all molecules in the system. It was used first in the mid-fifties of the last century [12] for two-dimensional modeling of hard disks system (2D-model of a monoatomic gas), and then was used tosimulate a variety of liquids, including water [13][14]. In 1972 this method was first applied to the simulation of a simplemodel of a linear polymer chain consisting of atoms connected by rigid bonds [15]. In 1975 the dynamics of short nalkanes was studied [16]. In subsequent years MD was used for detailed study of many specific molecules using both coarse-grained and detailed full-atomic models. The potential energy of these models usually include valence bonds, valence angles and dihedral angle energies as well as van der Waals and electrostatic energies. The definition of parameters set adequately describing the test molecule properties (force-field) is challenging and requires the experimental data for these molecules, quantum chemical calculations as well as iterative procedures and a very large amount of machine time. These calculations can be made only by large groups of specialists. Due to this reason several packages of standard computer programs, in which these parameters are defined for a fairly wide range of molecules become widely used in recent years. Currently the most popular molecular modeling packages are GROMACS, AMBER, CHARMM, and some others. Our simulation was performed by molecular dynamics method using the GROMACS 4.5.6 software package [17] and one of the most modern AMBER_99SB-ildn force fields [18].

Model and details of calculations
Modeling was performed using the molecular dynamics method for systems consisting of one lysine dendrigraft of second generation with positively charged NH3+ end groups, 16 molecules of KED peptide (with charge -1 each), water molecules and chlorine counterions in a cubic cell with periodic boundary conditions. The initial conformation for peptide with internal rotation angles of fi = -135º, psi = 135º, theta = 180º was modelled by Avogadro chemical editor. The structures were optimized in vacuum using molecular mechanics of AMBER force field. Further energy minimizations and simulations were performed using the GROMACS 4.5.6 software package and AMBER_99SB-ildn force fields. The potential energy of this force field consists of valence bonds and angles deformation energy, internal rotation angles, van der Waals and electrostatic interactions. The procedure of molecular dynamics simulation used for lysine dendrigraft and peptides has been described earlier (for dendrimers and linear polyelectrlytes) in [19][20][21][22][23][24][25][26][27][28][29]. In all calculations the normal conditions (temperature 300 K, pressure 1 ATM) were used. Computing resources on supercomputers "Lomonosov" were provided by supercomputer center of Moscow State University [30].

Results and Discussion
Snapshots of a system consisting of dendrigraft of second generation, peptides, ions and water during simulation are shown on Fig. 1 (water molecules are not shown for clarity). It is clearly seen that at the beginning of process (Fig. 1a) peptide molecules are far from dendrigraft. In the end of simulation (Fig. 1b) all peptide molecules are on dendrigraft surface. Atoms of dendrigraft molecule is shown as beads with diameter equal to their van der Waals radii. Valence bonds of various peptides are shown with lines of different colors (backbone of each peptide is shown by thick line of the same color as valence bonds). To characterize the size of the systems during the equilibration the mean squared gyration radius R g (t) was used. This function was calculated using g_gyrate function of GROMACS.

Complex formation
The time dependence of gyration radius R g at the beginning of calculation describes the process of compaction of subsystem consisting of dendrigraft and peptide molecules during formation of complex (see Fig.  2). From Fig. 2 it can be seen that complex of dendrigraft with 16 molecules of KED peptide forms within 20 ns. After this time the complexes size Rg fluctuate slightly, but its average value practically does not change with time. Therefore, we can assume that the systems become close to equilibrium state after 20ns.
Another quantity that can characterize the rate of complex formation is the total number of hydrogen bonds (N) between dendrigraft and peptide molecules. The dependence of this value on time is shown on Fig. 3 and demonstrates how the number of specific contacts between them increases during complex formation. This function was calculated using g_hbonds function of GROMACS.
From Fig. 3 it can be concluded that the number of Hbonds between dendrigraft and 16 molecules of KED peptide reaches equilibrium (plateau) for the first time after 20 ns. It correlates very well with the time dependence of the inertia radius shown in Fig. 2.  Fig. 3. Time dependence of hydrogen bond (Hbonds) number between DG2 dendrigraft and 16 molecules of KED peptide (b) during complex formation.

Equilibrium characteristics of complex
The sizes Rg of complex and dendrigraft in equilibrium state are evaluated by mean square of inertia radius averaged through the time t after equilibration (t>20ns). This values are given in table 1. In equilibrium state the sizes of the complex consisting of DG2 dendrigraft and 16 molecules of KED peptides (see Table 1) is greater than size of dendrigraft itself. It is quite natural, since it correlates with the molecular weight of the complexes increase compared to the molecular weight of the individual dendrigraft.
The shape of complex can be characterized by main components (Rg 11 , Rg 22 and Rg 33 ) of inertia tensor from Table 1 where Rg 11 is minimal and Rg 33 are largest and smallest eigenvalues of inertia tensor. The rough evaluation of anisotropy of dendrigraft and complex could be done using axial ratio Rg 33 / Rg 11 . This ratio is equal 1.43 for both dendrigraft and complex, It means that anisotropy of complex is practically the same as anisotropy of dendrigraft. It means that molecules of peptide are distributed evenly on dendrigraft surface in complex. Information about the internal structure of the equilibrium complex could be obtained using radial density distribution of atoms relatively center of inertia of system. These radial distribution functions (not normalized) are shown on Fig. 4. They were calculated using g_rdf function of GROMACS. Fig. 4 demonstrates that atoms of dendrigraft (curve 2) are located mainly in the center (close to distance r=0) of the complexes and atoms of peptides (curve 3) mainly on the surface of complex. They could only slightly penetrate into outer part of dendrigraft but not to its inner part (see curve 3 in Fig. 4). Distribution of all atoms belonging to complex (curve 1) is between distribution for dendrigraft atoms (curve2) and peptide atoms (curve 1).
The research is carried out using the equipment of the shared research facilities of HPC computing resources at Lomonosov Moscow State University [20]. I.M.N. and O.V.S. are thankful for the support from the grant of the Government of Russian Federation 08-08.