Controlling Peptide Function by Directed Assembly Formation: Mechanistic Insights Using Multiscale Modeling on an Antimicrobial Peptide–Drug–Membrane System

Owing to their potential applicability against multidrug-resistant bacteria, antimicrobial peptides (AMPs) or host defense peptides (HDPs) gain increased attention. Besides diverse immunomodulatory roles, their classical mechanism of action mostly involves membrane disruption of microbes. Notably, their unbalanced overexpression has also been associated with host cell cytotoxicity in various diseases. Relatedly, AMPs can be subject to aggregate formation, either via self-assembly or together with other compounds, which has demonstrated a modulation effect on their biological functions, thus highly relevant both for drug targeting projects and understanding their in vivo actions. However, the molecular aspects of the related assembly formation are not understood. Here, we focused in detail on an experimentally studied AMP–drug system, i.e., CM15–suramin, and performed all-atom and coarse-grain (CG) simulations. Results obtained for all systems were in close line with experimental observations and indicate that the CM15–suramin aggregation is an energetically favorable and dynamic process. In the presence of bilayers, the peptide–drug assembly formation was highly dependent on lipid composition, and peptide aggregates themselves were also capable of binding to the membranes. Interestingly, longer CG simulations with zwitterionic membranes indicated an intermediate state in the presence of both AMP–drug assemblies and monomeric peptides located on the membrane surface. In sharp contrast, larger AMP–drug aggregates could not be detected with a negatively charged membrane, rather the AMPs penetrated its surface in a monomeric form, in line with previous in vitro observations. Considering experimental and theoretical results, it is promoted that in biological systems, cationic AMPs may often form associates with anionic compounds in a reversible manner, resulting in lower bioactivity. This is only mildly affected by zwitterionic membranes; however, membranes with a negative charge strongly alter the energetic preference of AMP assemblies, resulting in the dissolution of the complexes into the membrane. The phenomenon observed here at a molecular level can be followed in several experimental systems studied recently, where peptides interact with food colors, drug molecules, or endogenous compounds, which strongly indicates that reversible associate formation is a general phenomenon for these complexes. These results are hoped to be exploited in novel therapeutic strategies aiming to use peptides as drug targets and control AMP bioactivity by directed assembly formation.


I. Extended Methods
I.1. All-atom simulation protocols Table S1. Applied restraints during the minimization and equilibration process of the all-atom simulations. Parametrization of Suramin. The suramin molecule was parametrized using the beta version of the MARTINI v3.0 1 force field. It was used instead of the latest MARTINI v2.2 2 due to its better performance during parameter validation when applied to the mapping.

Restraints
The all-atom reference data set was generated using the CHARMM36m 3 force field with the same simulation as for the other all-atom simulations. Suramin was simulated in the water phase with TIP3P water molecules, it was neutralized by Na + atoms, and the NaCl salt concentration was set to 150mM. The equilibration process was started by the steepest descent minimization in 5000 steps, which was followed by two restrained 25 ps long NVT simulations with Berendsen and velocity rescale coupling. Next, 450 ps NPT type equilibration was performed in three steps with Berendsen pressure coupling but different restraints on the suramin heavy atoms. As a final step of the equilibration, a 100 ns simulation was carried out with 200 kJ mol -1 nm -2 restraints on suramin heavy atoms using Nosé-Hoover thermostat with 1.0 ps coupling time and an isotropic Parrinello-Rahman barostat with 5 ps coupling to ensure the correct statistical ensemble. Finally, the 50 ns unrestrained production run was performed with the same parameters as in the last equilibration step. The step size was 2 fs, and the atomistic coordinates were written out in every ps, which resulted in 50000 frames for data fitting and validation. The coarse-grained representation of the molecule was created by an atom-to-bead mapping developed here (Figure 2, main text).
Based on the mapping scheme, the atom-to-bead mapped trajectory was created, and the initial bond, angle, and dihedral information was extracted and averaged over the whole simulation. The averaged values were then used as initial parameters for coarse-grained suramin. After building the initial structure, it was solvated with water molecules and neutralized by Na + ions where needed. The system was first simulated for one ns as an NVT ensemble (τt = 1.0 ps) coupled to a velocity rescaling thermostat with a target temperature of 300 K, which was followed by a 50 ns NPT simulation in 1 bar, controlled by an isotropic Berendsen barostat (τp = 6.0 ps, κp= 4.5 10 -4 bar -1 ). The trajectory was then analyzed and compared to the all-atom simulation results in terms of the bond, angle, and dihedral distributions, and an iterative process was applied to optimize the initial parameters using Python and R scripts.
Next, the production run was performed for 1 μs as an NPT ensemble, and velocity rescale thermostat was used for temperature coupling (τt = 1.0 ps) and Parrinello-Rahman barostat for pressure coupling (τp = 16.0 ps, κp= 3.0 10 -4 bar -1 ). The step size for every CG simulation was 20 fs.
The parameter validation process took place in two steps. First, the CG bond, angle, and dihedral parameters were compared directly, while in the second step, inherited descriptors such as the radius of gyration and octanol-water solvation free energy difference was matched with their all-atom counterparts.
The octanol-water solvation free energy difference was calculated using umbrella sampling 4 simulations and Plumed 2.4.3 software 5 . An octanol-water biphasic system was created where the suramin was in the center of the water phase neutralized by Na + atoms. The equilibration processes both for all-atom and CG simulations were the same as above. The collective variable of the umbrella sampling simulations was defined as follows: where is the z coordinate of the center of mass (COM) of the corresponding phase or molecule. Consequently, when = 1, suramin is in the center of the water phase, when = 0, it is inside the octanol phase. The following umbrella potential was used: where the force constant, = 5000 kJ/mol, and is the reference position of the window. The selected 0.04 step size resulted in altogether 26 umbrella windows. Each umbrella was simulated for 25 ns, the coordinates were saved every 25 steps, and the last 10 ns were analyzed. The Potential of Mean Force was calculated using the variational free energy profile (vFEP) approach 6 at 300 K temperature.
Validation of suramin parameters. The evaluation of suramin in a CG setup is described and briefly compared with its performance in related MD calculations, emphasizing its relationship with the phase surface present between a hydrophobic and a hydrophilic environment. The solvent phase allatom simulations were used as a basis for the suramin parametrization. As the mapped beads in the all-atom simulation defined narrow distributions in bond lengths and angles, the final bond, angle, and dihedral angle parameters were chosen accordingly. Consequently, most of the bonds were constrained, and stiff force constants were applied for the angles and for those bonds which were not constrained. Dihedral angles were defined to ensure the planarity of the sulfonated naphthyl rings, and the correct rotation of the two distinct parts of the suramin around the urea core. The distributions of the final parameters showed good overlap with their all-atom counterparts ( Figure S2-S25, end of this section).
Besides the direct comparison of the distributions of the parameters between the CHARMM36 allatom and the MARTINI v3.0b coarse-grained simulations, the parameters were also validated indirectly through two characteristics, the radius of gyration and the free energy difference of the water-to-octanol phase transition ( Figure S1). There is a significant overlap between the two distributions in terms of the radius of gyration; however, the expected value is lower in the coarse-grained simulation, which suggests a slightly more compact representation of suramin ( Figure S1A). While the radius of gyration is a good structural descriptor, it cannot measure the correctness of the chosen bead types; hence we calculated the free energy difference of the water-to-octanol phase transition, which reflects better on the merit of the electrostatic and van der Waals interactions. We chose the water-octanol system as it is commonly used for parametrization by calculating partitioning free energies 2 . The PMFs show approximately the same free energy difference for the transition in both cases, with one remarkable exception ( Figure S1B). Namely, the PMF of the transition in the CG simulation is a monotonic function, while it has a minimum at the phase boundary when simulated on an all-atom level. This indicates the "amphipathic" nature of suramin, which the coarse-grained parametrization could not capture. Nevertheless, the coarse-grained simulations show reasonable correspondence with their all-atom counterparts allowing its use in further simulations.