Interest to macromolecules with regular branched structure grows every year  . Beside well known dendrimers and in particular lysine dendrimers  and dendritic brushes  recently the new dendritic structure: dendrigraft was synthesized  . 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 consist entirely of lysine aminoacid residues   . At the same time their terminal groups could be functionalized by other aminoacids or by other active groups or molecules  . Lysine dendrigrafts are polymers that are rich with amines. Due to this reason they could be used as antibacterial  or antiviral agents  . It also makes possible the creation of biocompatible and biodegradable complexes with different drugs and other biological or bioactive molecules. In particular they could make complexes with oppositely charged peptides due to strong electrostatic interaction between their positively charged groups () and negatively charged amino acid 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 this ability to make complexes the dendrigrafts like other dendritic molecules could be used as multifunctional nanocarriers for delivery of drug or/and other molecules  for treatment of various disease.
Therapeutic Semax peptide  was selected for our study as a model peptide because it belongs to a class of regulatory peptides and has an antioxidant, antihypoxic and neuroprotective properties. Semax peptide is used for acute ischemic stroke prevention, during traumatic brain injury treatment, recovery of a patient after a stroke, in the case of optic nerve disease and glaucoma optic neuropathy. The drug is used in the form of solution for injection and as a spray. This peptide has molecular weight 863 Da and isoelectric point (pI) 5,13. Amino acid sequence of the peptide is shown on Figure 1.
In our previous paper we studied complexes of Semax peptides with lysine dendrimer  . The goal of this paper is to study the interaction between lysine dendrigraft of second generation and therapeutic peptide Semax using molecular dynamics method to determine whether the dendrigraft could adsorb peptide molecules and thus could be used for delivery of these peptides into cells.
2. Methods and Materials
2.1. Molecular Dynamics Method
Molecular dynamics (MD) method is currently the main method for simulation of
Figure 1. Amino acid sequence of Semax peptide.
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  for two-dimensional modeling of hard disks system (2D-model of a monoatomic gas), and then was used to simulate a variety of liquids, including water   . In 1972 this method was first applied to the simulation of a simple model of a linear polymer chain consisting of atoms connected by rigid bonds  . In 1975 the dynamics of short n-alkanes was studied  . In subsequent years MD was used for 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. 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  and one of the most modern AMBER_99SB-ildn force fields  .
2.2. Model and Calculation Method
Modeling was performed using the molecular dynamics method for systems consisting of one lysine dendrigraft of second generation containing 8 lysine residues in main chain, 48 positively charged groups in 8 side chains, 8 Semax peptides (with charge −1 each), water molecules and 40 chlorine counterions in a cubic cell (the edge of the cube equal to 9 nanometers) with periodic boundary conditions. The initial conformations for dendrigraft and peptides with internal rotation angles of j = −135˚, y = 135˚, q = 180˚ was constructed using molecular editor Avogadro. 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  -  . 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  .
3. Results and Discussion
Snapshots of a system consisting of dendrigraft of second generation, peptides, ions and water during simulation are shown on Figure 2 (water molecules are not shown Fclarity). It is clearly seen that at the beginning of process (Figure 2(a)) peptide molecules
(a) (b) (c)
Figure 2. Stages of the dendrimer and Semax peptides complex formation (initial, intermediate and final): system of dendrimer and 8 peptides: t = 0 (a); t = 30 ns (b); t = 200 ns (c). 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).
are far from dendrigraft. After 30 ns (Figure 2(b)) most of peptide molecules are already adsorbed on the surface of dendrigraft, and in the end (Figure 2(c)) all peptide molecules in both systems are on its surface.
To characterize the size of the subsystem containing dendrigraft and peptides during the equilibration the square of instant gyration radius was used. It was calculated using (1):
where R―is the center mass of dendrigraft, ri и mi―coordinates and masses of i atom correspondingly, N―is the total number of atoms in dendrigraft, M is the total mass of dendrigraft. This function was calculated using g_gyrate function of GROMACS software  .
3.1. Modeling of Equilibrium Process Establishment
The time dependence of radius of gyration Rg of dendrigraft and peptides at the beginning of calculation describes the kinetics of process of complex formation (Figure 3). From Figure 3 it can be seen that dendrigraft complex with 8 peptides forms within first 20 ns. After that the complex size Rg fluctuate slightly, but its average value practically does not change with time. Therefore, we can assume that the system is in equilibrium state.
The time dependence function of distances between dendrigraft and peptides also demonstrates the formation of complex within first 20 - 30 ns of simulation. Plateau on curves of Figure 3(b) (at times greater than 20 - 30 ns) means that system is in equilibrium state and all peptides are already adsorbed on surface of dendrigraft molecule. This function was calculated using g_bond function of GROMACS software  . Thus the dendrigraft makes complex with 8 peptides within first 20 - 30 ns of simulation.
Another quantity that can characterize the rate of complex formation is the total number of hydrogen bonds (N) between dendrigraft and peptides. The dependence of this value on time is shows on Figure 4 and demonstrates how the number of specific contacts between them increases during complex formation. This function was calculated using g_hbonds package of GROMACS.
Figure 3. System of dendrigraft DG2 and 8 Semax peptides: (a) time dependence of radius of gyration, (b) time dependence of distance between dendrigraft and peptides.
Figure 4. Time dependence of hydrogen bonds number (N) during the complex formation: DG2 and 8 Semax.
3.2. Modeling of the Equilibrium State
The size of complex in equilibrium state is evaluated by mean square of radius of gyration averaged through the time after equilibration (2):
where Dt = tmax − teq , tmax―full time of calculation and and teq―time of equilibration and < > means averaging the through equilibrium part of MD trajectory.
In equilibrium state the size of the complex (DG2 and 8 Semax peptides) is larger than the size of dendrigraft (see Table 1). It is quite natural, since it correlates with the increase of molecular weight of the complexes compared to the molecular weight of the individual dendrigraft.
The shape of complex can be characterized by its tensor of inertia main component ratio (, ,), that are in Table 1. For example, in the simplest case, anisotropy can be characterized by.
The shape of complex could be roughly characterised by ratio of largest and smallest eigenvalues of inertia tensor describing our system. Calculated values of these anisotropy for our systems are presented in Table 2. Thus the anisotropy of complex in water is slightly less than anisotropy of single dendigraft in water itself.
Information about the internal structure of the equilibrium complex could be obtained using radial density distribution (3) of different groups of atoms relatively center of inertia of system.
where mcomp―mass of all atoms in complex; Vкомп―volume of complex.
Table 1. Eigenvalues, , of tensor of inertia in dendrigraft without or with peptides.
Table 2. The values of anisotropy of shape for dendrigraft and for complex with peptides.
Figure 5. Radial distribution p(r) curves: dendrigraft DG2 and 8 Semax. Distribution curves: peptide atoms (1); dendrigraft atoms (2); all atoms of complex (3).
Figure 5) is located mainly in the center of the complexes and peptides (curve 1, Figure 5) mainly on the surface of complex. At the same time some fraction of peptides could slightly penetrate into outer part of dendrigraft but not to its inner part (see Figure 5).
The number of hydrogen bonds between peptides and dendrigraft shows how tightly peptides associate with dendrigraft. From equilibrium part of Figure 4 (at times greater 20 - 30 ns) it follows that average hydrogen bonds number in equilibrium state between dendrigraft DG2 and 8 Semax peptides is equal to 25.
The distribution function of hydrogen bonds number (Figure 6) shows how the number of hydrogen bonds in the equilibrium state can deviate relative to the average value. We obtained that the resulting function for complex has a peak of numbers of bonds equal 25 which coincide with average number of hydgrogen bonds between dendrimer and peptides. Thus this function is quite symmetrical. Fluctuations in hydrogen bonds number for our system are in the range of 15 - 40.
The other characteristic of interaction between dendrigraft and peptides in complex is the distribution of ion pairs number between their charged groups. Figure 7 shows the dependence of probability (not normalized) of different distances between positive groups of dendrimer and negative groups of peptides. It is seen that there is a sharp first peak, corresponding to the direct contact between positively charged groups () of dendrigraft and negatively charged groups (COO−) of the glutamic acid in peptides. There is also small second peak on this dependence probably due to second neighboring charges.
Figure 6. The distribution function P(N) of hydrogen bonds number N of complex: complex DG2 and 8 Semax.
Figure 7. Radial distribution function of ion pairs:complex of DG2 and 8 Semax molecules, groups of dendrimer and COO− groups of peptides (1); groups of dendrimer and Cl− ions (2).
To evaluate the translational mobility of our systems, the time dependence of the mean square displacements (4) of the centers of inertia (MSD), were calculated (Figure 8). MSD was calculated using g_msd function of GROMACS  .
We have found that dependence of MSD function on time is almost linear in some interval of time t in double logarithm coordinates (not shown). It means that in this interval the motion of complex is the diffusion-like motion (see Figure 8). Coefficient of translational diffusion of the complex was determined from the slope of the obtained time dependence of MSD for system and was equal D = 0.18 ± 0.01 sm2/s.
The system consisting of second generation dendrigraft and 8 Semax molecules in water with counterions has been studied. The process of complex formation by lysine dendrigraft of second generation and therapeutic Semax peptides and the equilibrium structure of the complex was investigated for the first time by the method of molecular dynamics simulation. It was shown that dendrigraft-peptide complex formation occurs rather quickly (in first 20 - 30 ns). After this time the complex is stable and the dendrigraft atoms are mainly inside the complex, while the peptide atoms are mainly on the surface of dendrigraft. In future it is necessary to study process of complexes formation and stability of similar complexes with greater number of peptides to define maximum capacity of this particular dendrigraft for delivery of Semax peptide.
This work was partly supported by grant 074-U01 of Government of RF and RFBR
Figure 8. Mean square displacements of the centers of inertia: complex of DG2 and 8 Semax.
grants 16-03-00775 and 15-33-20693mol_a_ved. Computing resources on supercomputers “Lomonosov” were provided by supercomputer center of Moscow State University  .