Uncovering the mechanisms of cyclic peptide self-assembly in membranes with the chirality-aware MA(R/S)TINI forcefield

Uncovering the mechanisms of cyclic peptide self-assembly in membranes with the chirality-aware MA(R/S)TINI forcefield Alfonso Cabezón a , Martin Calvelo a,b , Juan R. Granja a , Ángel Piñeiro c, ⇑ , Rebeca Garcia-Fandino a, ⇑ a Department of Organic Chemistry, Centro Singular de Investigación en Química Biolóxica e Materiais Moleculares (CiQUS), Santiago de Compostela University, CIQUS, Spain b Departament de Química Inorgánica i Orgànica and Institut de Química Teòrica i Computacional (IQTCUB), Universitat de Barcelona, Barcelona 08028, Spain c Soft Matter & Molecular Biophysics Group, Department of Applied Physics, Faculty of Physics, University of Santiago de Compostela, Spain graphical abstract article info Article history: Received 2 February 2023 Revised 19 February 2023 Accepted 16 March 2023 Available online 22 March 2023 Keywords: Cyclic peptides Coarse-grained molecular dynamics simulations Martini forcefield Lipid bilayers abstract Cyclic peptides (CPs) formed by alternation of Dand L-amino acids (D,L-CPs) can self-assemble into nanotubes (SCPNs) by parallel or/and antiparallel stacking. Different applications have been attributed to these nanotubes, including the disruption of lipid bilayers of specific compositions and the selective transport of ions throughout membranes. Molecular dynamics (MD) simulations have significantly contributed to understand the interaction between CPs, including the structural, dynamic and transport properties of their supramolecular aggregates. The high computational cost of atomic resolution forcefields makes them impractical for simulating the self-assembly of macromolecules, so coarse-grained (CG) models might represent a more feasible solution for this purpose. However, general CG models used for the simulation of biomolecules such as the MARTINI forcefield do not explicitly consider the noncovalent interactions leading to the formation of secondary structure patterns in proteins. This becomes particularly important in the case of CPs due to the Dand L-chirality alternation in their sequence, leading to opposite orientations of the backbone polar groups on both sides of the cyclic ring plane. In order to overcome this limitation, we have extended the MARTINI forcefield to introduce chirality in each residue of the CPs. The new parametrization, which we have called MA(R/S)TINI, reproduces the expected selfassembly patterns for several CP sequences in the presence of different membrane models, explicitly considering the chirality of the CPs and with no significant extra computational cost. Our simulations provide new mechanistic information of how these systems self-assemble in presence of different lipid scenarios, showing that the CP-CP and CP-membrane interactions are sensitive to the peptide sequence chirality. This opens the door to design new bioactive CPs based on CG-MD simulations. A web-based tool for https://doi.org/10.1016/j.jcis.2023.03.101 0021-9797/Ó2023 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ⇑ Corresponding authors. E-mail addresses: [email protected] (Á. Piñeiro), [email protected] (R. Garcia-Fandino). Journal of Colloid and Interface Science 642 (2023) 84–99 Contents lists available at ScienceDirect Journal of Colloid and Interface Science journal homepage: www.elsevier.com/locate/jcis the automatic parameterization of new CP sequences using MA(R/S)TINI, among other functionalities, is under construction (see http://cyclopep.com). Ó2023 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction Molecular self-assembly is the spontaneous organization process of randomly distributed and oriented molecules into welldefined and stable structural arrangements without guidance from an artificial external source [1]. This is mainly controlled by the cooperative result of many weak non-covalent interactions of different type: ranging from electrostatic, direct or water-mediated hydrogen bonds (H-bonds), halogen bonds and van der Waals interactions among others. Nature, through billions of years of evolution, is an excellent reference to find inspiration and ideas for molecular assembly [2]. Unsurprisingly, many scientists have found motivation in the captivating and sophisticated biological recognition mechanisms and they are making an effort to mimic them with diverse organic and inorganic building blocks in the laboratory, with the aim of creating new materials, molecular carriers for active principles or specific molecular traps for toxic compounds, amongst other applications. Peptides represent one of the most evident organic building blocks for the design of nanomaterials with controllable structural features, due to their simple topology, relative chemical and physical stability, diversity in sequences and shapes, biocompatibility, biodegradability and feasibility to synthesize in large amounts [3,4]. In particular, cyclic peptides (CPs) are a fascinating class of molecules that can be designed, because of their conformationally restricted freedom, to fold or self-assemble into diverse monoand multidimensional suprastructures with potential applications in biomedicine, nanoelectronics, or catalysis [5–7]. As bioinspired supramolecular building blocks, CPs can stack into selfassembled cyclic peptide nanotubes (SCPNs), which have attracted special attention due to the facility to modify their structure and, thus, their functional properties [5,7]. This makes them well suited to be used in numerous applications, such as electronic devices, artificial photosystems, photoresponsive materials, biosensors, antimicrobials and antiviral agents, selective transmembrane transport channels, catalysis, and drug delivery. The first SCPNs were proposed in 1974 [8], but they could not be synthesized until 1993 [9]. These materials comprised the use of D,La -cyclic peptides (D,La -CPs) formed by an even number of alternating Dand La -amino acids. This leads to a flat ring-shaped conformation in which amide groups (NH and CO) lay perpendicular to the plane of the CPs, allowing the formation of supramolecular nanotubes upon hydrogen bonding with neighbor CPs in a b-sheet fashion (parallel or antiparallel, Fig. 1). Such a special arrangement leaves an empty inner duct of the assembly, with all the side chains exposed on the external surface of the cylindrical structure. One of the main advantages of these nanostructures is the simplicity with which their internal and external properties can be modified based only on the selection of the appropriate CP sequence: the inner diameter of the SCPNs is determined by the number of residues, whereas their outer surface properties depend on the amino acid side chains. The properties and orientation of the SCPN can be tuned by tailoring the peptide sequence [10]. In this sense, SCPNs with appropriate hydrophobic sequences have emerged as attractive transmembrane channel mimics that are able to replicate specific functions of natural transport systems, in terms of affinity, efficiency, stability, and selectivity [11,12]. Alternatively, amphipathic D,La -CPs form SCPNs that lie parallel to lipid membranes causing their disruption, thus making these materials particularly suitable for the development of new antimicrobial agents (Fig. 1A) [13–15]. Molecular self-assembly and, in particular, the assembly of CPs into SCPNs, typically involves a complex interplay of physical interactions that is extremely challenging to investigate with traditional experimental techniques. Different simulation methods can be employed nowadays to model CPs and SCPNs thanks to the currently available computational facilities [16,17]. In particular, Molecular Dynamics (MD) simulations allow to reproduce the self-assembly process of this kind of molecular systems, so it can facilitate unprecedented access to a range of physical properties at the intra-, inter-, and supramolecular scales. While atomic resolution models (AT) can provide important insights into SCPNs [18– 25], they are too computationally expensive to follow dynamic behavior of the self-assembly process of CPs. To tackle these limitations, coarse-grained (CG) models, in which multiple atoms are grouped to form joint interaction sites, can be used. Thus, some intramolecular degrees of freedom are ignored to increase the sampling in the system dimensions leading to the spontaneous formation of supramolecular structures. Another advantage of using CG resolution is that the time step employed to solve the motion equations in MD simulations can be significantly increased without losing stability in the calculations. This is because the fastest intramolecular movements are no longer considered at this resolution level, and the shortest periods of the remaining degrees of freedom are typically longer than 100 fs. The immediate consequence of using longer time steps is that significantly longer trajectories, compared to those obtained using AT resolution, can be generated with the same computational resources. Multiscale approaches to take the best from both resolution levels -CG for self-assembly and AT to identify atomic level interactions stabilizing the supramolecular patternshave been proposed [26]. MARTINI is the most popular CG forcefield for the simulation of biomolecular systems [27–29]. It has been used to investigate the spontaneous self-assembly of peptides, proteins, carbohydrates, surfactants, and polymers into different supramolecular structures. Despite its unquestionable advantages, i.e., reduced computational costs and accurate description of molecular movements, none of the available versions of this forcefield [30–33] considers the non-covalent interactions leading to the formation of secondary structures in peptides and proteins. Secondary structure arrangements are imposed by harmonic potentials between backbone beads. Thus, MD simulations based on MARTINI are not suitable for the study of transitions among different secondary structures, like for example in protein folding. A direct consequence of this limitation is the loss of chirality in the amino acids, and also the lack of information about directional non-covalent forces such as hydrogen, halogen or chalcogen bonds [34], which becomes particularly important in the case of CPs due to the Dand L -chirality sequential fluctuation: the donor and acceptor groups of each residue (NH and CO) are alternatively oriented in opposite directions with all amide groups of the residues of identical chirality pointing to the same side of the CP plane Due to this directional orientation, D residues of one CP unit interact with L residues of a contiguous, parallel CP, while residues of the same chirality interact in contiguous, antiparallel CPs. (Fig. 1A). Consequently, the self-assembly of D,La -CPs using MARTINI produces a misleading representation of reality because it does not consider these specific interactions, thus making it impossible to distinguish between the parallel A. Cabezón, M. Calvelo, J.R. Granja et al. Journal of Colloid and Interface Science 642 (2023) 84–99 85 and antiparallel b-sheets that can be formed along the SCPN (Fig. 2) [7]. Furthermore, the number of possible interactions between backbone beads of different CP units is overestimated, allowing stackings that should be forbidden due to the previously mentioned alternated chirality and alignment (Fig. 2). This overestimation of these interactions is extended not only to all the b-sheet based protein/peptide structures, but also to all the systems where the directionality of the H-bonds is critical, even if the residues have the same chirality. All the CG-MD simulations that have been carried out so far employing MARTINI are subjected to these limitations [13–15,35–37]. In order to solve this problem, an extension of the MARTINI forcefield that introduces chirality in the residues is proposed here. While ad hoc CG models might lead, in principle, higher accuracy in Fig. 1. A. Top-left: general structure of a CP, with an even number of alternating D- (in blue) and L- (in red) a -amino acids. Top-right: structure of parallel and antiparallel dimers formed by two CPs, as models of the basic interactions between stacked CPs together with a schematic representation of a SCPN (stacked yellow rings). Bottom: schematic representation of a transmembrane channel formed by the self-assembly of hydrophobic CPs (left) and of a SCPN parallel to the membrane formed by amphipathic CPs (right). B.Top: Structure of the three simulated CPs, CP1,CP2 and CP3. The underlined residues in the abovementioned sequences represent the non-natural D -residues. Bottom: CG representations of the membrane models M1 (mammalian), M2 and M3 (bacteria), used in this work. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) A. Cabezón, M. Calvelo, J.R. Granja et al. Journal of Colloid and Interface Science 642 (2023) 84–99 86 the representation and treatment of a specific system, we preferred to extend a well stablished model for which many functional groups have already been parameterized. The new parametrization, which we have called MA(R/S)TINI (Fig. 2), reproduces the expected self-assembly patterns for several CP sequences, explicitly considering the chirality of the CPs with no significant extra computational cost. This strategy facilitates monomer customization as well as the simulation of CPs with relatively complex molecular mixtures: lipid bilayers of different composition, proteins, linear peptides, surfactants, etc. In this way, the selfassembling of monomeric peptides with different sequences can be studied in a variety of environments, allowing to consider structure-activity relationships. The self-assembly process of three different CP sequences in the presence of membrane models for healthy mammal and for bacteria has been simulated using the classical and the new parameterizations of the MARTINI forcefield (Fig. 1B). The three CPs are expected to exhibit different behaviors in the presence of the employed membrane models due to the large differences in the total charge, as well as in the hydrophobic and electrostatic dipolar moments. The MA(R/S)TINI parameterization reproduces the general trends expected for the proposed sequences while it explicitly considers the chirality of the CPs that is missed in the original MARTINI model. Furthermore, the simulations with the new parameterization provide new mechanistic information of how these systems self-assemble in presence of different membranes, at a level of resolution that had not been previously attained. 2. Materials and methods 2.1. Peptide sequences and membrane models Three cyclic sequences of 8 amino acids each, with different expected behaviors in the presence of a lipid environment, were used (Fig. 1B): CP1 (RRKWLWLW), with amphipathic character and with demonstrated activity against bacteria [38];CP2 (LLLLLLLL), consisting on a fully hydrophobic sequence so, a priori, it should lead to the formation of transmembrane SCPNs; and CP3 (LQLWLWLW), mostly hydrophobic and known to form transmembrane nanotubes [11]. The underlined residues in the abovementioned sequences represent the non-natural D -residues. These CPs have been simulated in the presence of three different lipid membrane models, depicting a mammalian and two bacteria cell membranes (Fig. 1B). The mammalian model (M1) consisted of a pure 1,2-dimyristoil-sn-glycero-3phosphatidylcholine (DMPC) bilayer, as this lipid represents >50% of the phospholipids in most eukaryotic membranes [39]. The bacterial membrane models (M2 and M3) were formed by a mixture of 1,2-dimyristoyl-sn-glycero-3-phospho-(1 0 -rac-glycerol) (DMPG) and DMPC in 3:1 ratio for M2 and 1,2-dimyristoyl-sn-glycero-3phosphoethanolamine (DMPE) and DMPG in 1:9 ratio for M3, since these lipids are the main constituents of Gram-Negative and GramPositive bacteria membranes [40,41]. These models are expected to provide useful information about the role of lipid composition, and Fig. 2. Atomistic representation of a generalized CP backbone with the optimal antiparallel and parallel conformations of CP dimers, along with a representation of the backbone in the MA(R/S)TINI extension of the force field and the standard MARTINI. Note that standard MARTINI does not account for the limitation of certain orientations, making it impossible to distinguish between parallel and antiparallel b-sheets that can be formed along the SCPN. The ‘‘Unmatched H-bonds!” labels indicate that with those relative rotations of both CP units (45°) the formation of H-bonds is not feasible. A. Cabezón, M. Calvelo, J.R. Granja et al. Journal of Colloid and Interface Science 642 (2023) 84–99 87 in particular of the net charge of the membrane, in the selfassembly of the CPs. 2.2. Structure and parameterization of the CPs The coordinates of the particles describing the CP models and the corresponding force field parameters (topologies) were generated as follows: 2.2.1. Standard MARTINI Atomistic coordinates of linear peptides with the required sequences are automatically obtained using the Python PeptideBuilder library [42]. The resulting structures are then mapped into a MARTINI v2.2 CG model [43] with the martinize.py tool [31]. Next, following the Tarek’s strategy [37], the topology provided by this tool is modified by introducing the bonds, angles and dihedrals that reproduce a reliable structure for the CPs: The bonds between backbone beads (BB) are described by a harmonic potential with the equilibrium distance at 3.8 Å and a force constant of 6275 kJmol 1 nm 2 . The angle between consecutive BB beads is described using again a harmonic potential with an equilibrium angle h 0 of 135°and a force constant of 627 kJmol 1 rad 2 . Improper dihedrals centered on each BB and using the corresponding lateral chain together with the two bound BB beads are imposed to restrain the relative orientation of each residue. The angle and force constants employed in this case are n 0 = 180° and K n = 418 kJmol 1 rad 2 . A second group of improper dihedrals between 4 consecutive BB beads with an equilibrium angle of 0° and the same force constant employed above (418 kJmol 1 rad 2 ) is added to restrain the planarity of the CPs. All these parameters were directly taken from the Tarek’s previous work [37], with the exception of the force constant for the bonds between BB beads, In order to increase the rigidity of these interactions and enhance the stability of the simulations, such force constant was augmented by a factor of 100. The resulting structures and parameters were taken as a reference to include the chirality in the extended version of the forcefield (see below). 2.2.2. MA(R/S)TINI parameterization To introduce chirality in the Dand Lresidues composing the CPs, the original P5 bead representing the backbone in the classical MARTINI CG model was replaced by a new particle (CPBB) combined with two virtual sites (CPCO and CPNH). The virtual sites (VS) were introduced to depict the C@O acceptor (CPCO) and NAH donor groups (CPNH) of the amino acids (Fig. S1). The parameters associated to the new beads ( r and e values for the LennardJones potential, equilibrium distances, angles, and force constants) were optimized through an extensive trial and error iterative process in which several tens of combinations were essayed. The criterion utilized for the determination of the optimal LJ parameters was guided by a purely empirical approach, which involved screening a large number of potential combinations until the most favorable ones were identified. The final values lead to both parallel and antiparallel SCPNs with the experimentally reported distances between CPs and the expected interactions between the VS (Fig. 2)[9,11]. This extension of the MARTINI force field was designed to maintain the versatility of its original version, so the modifications only affect the interaction between CPs. Thus, the new CPBB bead has exactly the same parameters as the original P5 MARTINI particle, including the cross-interactions with all the beads of the force field. However, the r Lennard-Jones (LJ) parameter for the CPBB-CPBB self-interaction was decreased from 0.47 to 0.4 Å. This reduction was key to get a CP-CP distance closer to the experimentally reported by electron diffraction and Fouriertransformed infrared spectroscopy values [9,11], compared to the original version of the force field (see below). On the other hand, VS were used to model the C@O and NAH groups, thus introducing chirality in the CG model of the CP, directly keeps them at fixed positions with respect to the backbone. 3out type VS, as defined in GROMACS [44], were employed for both CPCO and CPNH. This means that the position of each VS is given by a set of coordinates in a local reference system given by three beads. The position in the plane and out of the plane defined by such beads can be easily set in the topology file. Trying to mimic the optimized atomic structure of the CPs, both VS particles, CPCO and CPNH, were initially placed at 1/3 of the distance between consecutive BB beads (d), along the axis joining them (Fig. S1). This distance was set to 0.35 dupon optimization. A similar process was followed to get the distance to the central bead (CPBB) in the vertical axis (h), obtaining a value of 1.053 Å (which is the length of the NAH bond in atomic resolution models) for both VS. This distance, combined with the reduced r value of the CPBB beads, is key to define the distance between CPs in SCPNs. In contrast to CPBB, the new CPCO and CPNH particles are allowed to interact only between them. The self-interaction parameters are r = 0.35 nm and e =3kJmol 1 , leading to a relatively weak interaction, while the crossinteraction is much stronger: r = 0.26 nm and an e = 5.6 kJmol 1 . 2.3. Construction of the lipid bilayers CG lipids distributed in 250 units per leaflet were introduced in a 12.5x12.5x15 nm 3 rectangular simulation box and solvated using the MARTINI Builder module of the CHARMM graphical-user interface [45] (Fig. 1B). 2.4. Free MD simulations Once the CP structures were optimized for both parametrizations (standard MARTINI v2.2 and MA(R/S)TINI), different number of CPs (20, 40, 80 and 160) were inserted into the water phase of the simulation boxes corresponding to each of the three membrane models. As a result, 36 systems (3 membrane models 3CP sequences 4 CP concentrations) were prepared. The counterions required to neutralize the electrostatic charge of the whole system were added in each case. A steepest descent energy minimization was performed to remove eventual atomic overlaps and clashes, followed by 50 ps (for simulations in MARTINI) and 50, 200, 2500, and 1500 ps (for MA(R/S)TINI) of stepwise constant pressure equilibrations using the Berendsen barostat [46]. The production trajectories were carried out at 1 bar and 303 K using the semiisotropic Parrinello–Rahman barostat [47], and a V-rescale thermostat [48]. The LINCS algorithm [49] was used to eliminate bond vibrations or its drifting in the system. The Particle Mesh Ewald [50] method with periodic boundary conditions was used to treat the long-range electrostatics with a 1.1 nm cut-off in direct-space. The dielectric constant was set to 15, as recommended for the non-polarizable version of MARTINI [51]. Van der Waals interactions were computed using a spherical cut-off of 1.1 nm. Production trajectories of 20 l s with a time step of 25 fs (for the MARTINI model) and 20 fs (for the MA(R/S)TINI model) for the integration of the motion equations were performed. Energies and positions of atoms were stored every 1000 steps (25 ps). All MD simulations were performed with the GROMACS 2021.4 [52] software package. 2.5. Metadynamics simulations Parallel and antiparallel CP dimers were built for both MARTINI parameterizations (standard v2.2 and MA(R/S)TINI versions of the forcefield). Each system was subjected to a 500 ns free MD simulation following the protocol described above, in all cases the dimers remain stable. Then, well-tempered two-dimensional metadynamA. Cabezón, M. Calvelo, J.R. Granja et al. Journal of Colloid and Interface Science 642 (2023) 84–99 88 ics simulations with PLUMED 2.8.0 [53] and GROMACS 2021.4 [52] were performed at 298 K, using the distance between the centre of mass of both dimers and the relative angle between the vectors joining the centre of each CP with its own first BB bead as collective variables. This angle measures the relative rotation between CPs while they share the same symmetry axis (perpendicular to the plane of the peptide ring). Gaussians with 0.1 kJ mol 1 height and 0.01 nm width were deposited every 500 ps using a biasing factor of 15. Several harmonic restrains were applied to the relative movement of both CPs: the angle between the normal to the CP rings was set to 0°for parallel dimers and to 180°for antiparallel structures; and the centres of mass were restrained to move throughout the axis perpendicular to both CPs, with no lateral displacement. The force constants used for both restraints was 1000 kJ mol 1 rad 2 for the angle and kJ mol 1 nm 2 for the distance to the axis. An artificial wall was also set to prevent a separation longer than 2.5 nm between the centre of both CPs. A minimum of 10l s-long trajectories with 6 walkers were obtained following this protocol. 2.6. Analysis of the trajectories VMD [54] was employed to generate snapshots and animations from the MD trajectories. The analysis of the simulations was carried out using GROMACS and PLUMED tools as well as specific code written in Python mainly based on the MDAnalysis [55], NumPy [56] and Matplotlib libraries [57]. For the free simulations, the number of CPs forming nanotubes of different size was computed as a function of time. For the kinetic analysis, two CPs are considered to form a dimer when the distance between their centres of mass is lower than 6 Å. This analysis allows to observe the spontaneous self-assembly of independent CP molecules, as well as their interaction with the different membrane models, since they are initially at random positions and orientations within the simulation box. Average values of the following properties along the last microsecond of the trajectories were determined: the number of CPs forming parallel and antiparallel dimer b-sheets, the distance between the centers of mass of the same CPs, and their rotational angles. For this analysis a more restrictive cutoff of 5.5 Å was employed to consider that two contiguous CPs form a dimer. The parallel and antiparallel orientations of contiguous molecules are identified from the angle between the normal vectors to the CP rings (Fig. S2). These normal vectors are determined considering the closest plane to the 8 CPBB beads of each CP, and they are clockwise oriented in the direction defined by the residue sequence. Then, the angle for parallel and antiparallel relative orientations are expected to be 0°and 180°, respectively. Potential of mean force (PMF) profiles were obtained from the 2D metadynamics simulations using the sum_hills tool of PLUMED [53]. 3. Results and discussion 3.1. Validation of the MA(R/S)TINI parameters: Potential of mean force profiles for parallel and antiparallel CP dimers In order to evaluate the reliability of the MA(R/S)TINI parameters, metadynamics simulations (see Methods) were carried out for parallel and antiparallel dimers with the three selected CPs (CP1-CP3). The objective of the metadynamics simulations is to quantify the energy of interaction between CP dimers as a function of the distance between the center of mass of CP units, of the rotational angles and also of the relative orientation between them (parallel or antiparallel). The MA(R/S)TINI parameters led to significantly more reliable results than the original MARTINI parameterization in all the studied cases. The resulting energy profiles (Fig. 3) clearly confirm: (i) the minimum energy for all peptides and for both the parallel and the antiparallel orientations corresponds to a shorter distance between CP units using MA(R/S)TINI than using MARTINI; (ii) only 4 rotational angles centered at 0°,90°, 180°and 270°are allowed for the new chirality aware parameterization while for the original MARTINI parameterization 8 angles are observed; (iii) the rotamers observed for the original force field are shifted by 22.5°with respect to the expected angles; and (iv) in contrast to the new parameterization, MARTINI is not sensitive to the relative orientation (parallel or antiparallel) of the CP units. In addition to this information, several interesting features can be seen in the PMF profiles. First, both parameterizations exhibit periodic energy barriers of decreasing amplitude as the distance between both CPs increases. Such energy barriers arise from the solvation layers with the explicit water beads considered in the force field, which are identical for both parameterizations (Fig. S3). This behaviour is reflected in the 2D PMF profiles as vertical bands of different intensity. In the 1D projection of the PMF on the distance between CPs, the solvation layers are seen as damped waves with one or more small shoulders along the first period (Fig. 3). Those shoulders are provoked by the impact of the rotamers in the first solvation layer. Such effect is clearer for the original MARTINI probably because for that parameterization there are twice more rotamers than for MA(R/S)TINI. This behaviour is clearly observed for the three CPs and for both relative orientations (Figs. S4-S8). Notably, MA(R/S)TINI distinguishes better the sequence than MARTINI. More specifically, the depth of the minima for CP1 and CP3 are more similar to each other than those for CP2. This happens for both parameterizations and relative orientations, but the differences are clearer for the new parameterization (see Table 1). The energy values in Table 1 also show that the parallel orientation is slightly preferred over the antiparallel one when using the MA(R/S)TINI model while the PMF profile is much more similar for both relative orientations of the CP units for old parameterization of the force field. The anisotropy of the rotations is only marginal in the metadynamics simulations in pure water, in contrast to what is observed in the unbiased free MD simulations (see section 3.2.3). This indicates that the external presence of the membrane model favours specific rotations of the CPs, more than the interaction between contiguous CPs. 3.2. Self-assembly of CPs into SCPNs Once it is confirmed that MA(R/S)TINI reproduces the correct relative rotations between CP dimers, unbiased MD simulations using 20, 40, 80 and 160 CP units in the presence of three different membrane models (M1-M3), were carried out for all the sequences (CP1-CP3). The results of the analysis obtained from all these simulations are presented in what follows. The presence of lipid membranes in the simulation boxes represents a double-edged sword for the self-assembly of CPs into SCPNs. On one hand, the lipid membranes may compete with CP aggregates for the capture of peptide monomers, hindering the formation of SCPNs in aqueous media. On the other hand, the lipid bilayers can also help the self-assembly of CPs by attaching individual molecules to their surface, increasing their local concentrations and potentially aiding the formation of SCPNs embedded in membranes. The interaction mechanism of the CPs or the corresponding nanotubes with the lipids will be determined by the nature of the amino acids in each sequence. As mentioned in the Methods section, CP1 is expected to form horizontally oriented SCPNs in the lipid membrane due to its amphipathic character, while CP2 and CP3, are expected to form transmembrane nanotubes because of their hydrophobic sequences. A. Cabezón, M. Calvelo, J.R. Granja et al. Journal of Colloid and Interface Science 642 (2023) 84–99 89 The membrane models, CP structures and concentrations chosen for the free simulations are intended to provide insight into the ability of both MARTINI and MA(R/S)TINI force fields to accurately represent not only the ability to form nanotubes but also to anticipate trends for each system. The kinetic variation in cluster sizes over the whole trajectories, as well as the distributions Fig. 3. Two-dimensional potential of mean force (PMF) profiles obtained for the interaction between CP2 dimers with antiparallel relative orientation from metadynamics simulations (see Methods section) using the MA(R/S)TINI (top) and the MARTINI (bottom) models. The 1D projections of the PMF on each collective variable (distance between CP units and rotational angle) calculated over periods of 0.8 l s is shown on top and on the right of the PMF plots. The results obtained for all the sequences and orientations are shown in Figs. S4-S8. Table 1 Energies (in kJ/mol) and their corresponding standard deviation (between parenthesis), corresponding to the flat region in the one-dimensional projection of the PMF profile on the first collective variable (the distance between monomers) for the 3 CPs with both relative orientations and parameterizations of the force field. These energies represent a measurement to the depth of the lowest minima with positive sign. MA(R/S)TINI MARTINI Parallel Antiparallel Parallel Antiparallel CP1 159 (17) 115.6 (7.0) 76.7 (5.3) 75.7 (5.4) CP2 115 (16) 101.9 (7.0) 67.1 (4.0) 66.7 (3.4) CP3 154 (21) 114.6 (7.6) 78.5 (5.0) 78.6 (5.9) A. Cabezón, M. Calvelo, J.R. Granja et al. Journal of Colloid and Interface Science 642 (2023) 84–99 90 between parallel and antiparallel CP arrangements, CP-CP distances and rotational angles for contiguous CPs in the last microsecond of each trajectory, will be used as quantitative descriptors of the self-assembly process and the resulting final supramolecular structures. 3.2.1. Kinetics of the simulations and SCPN length distributions In all cases, the distribution of CPs in aggregates of different sizes was stable before the end of the trajectories (20 l s long). This analysis revealed significant differences based on the peptide sequence and on the lipidic environment sensed by the CPs (Fig. 4 and S9-S15). The results of the simulations in the presence of different membrane compositions are discussed separately in the next section: 3.2.1.1. Simulations in the presence of DMPC bilayers (M1). The results of the simulations containing 40 CP units and using the MA(R/S)TINI model are shown in Fig. 4. It is observed that, even after 20 l s, a significant number of CP1 and CP3 monomers remain as independent entities in the corresponding simulation boxes, while for the simulation with CP2 almost all the peptides become part of trimers or larger aggregates within the first 7 l s. The distribution of CP1 molecules between SCPN of 1–5 molecules (monomers to pentamers) is quite homogeneous, ranging from 5 to 10 units per type of assembly over the last 10 l s of the trajectory. The last few l s contain 6 monomers, 5 dimers, 2 trimers, 2 tetramers and 2 pentamers. CP2 clearly forms larger and very stable aggregates, the last 14 l s with 2 dimers, 5 trimers, 1 tetramer, 2 pentamers and even 1 heptamer, with no monomers at all. The behaviour of CP3 is more similar to CP1 than to CP2, the last few l s containing 9 monomers, 8 dimers, 2 trimers, 1 tetramer and 1 pentamer. Similar general trends, although with some singularities, are observed at other peptide concentrations. CP2 still exhibits a higher ability to aggregate, compared to CP1 and CP3, in all the simulations. For the simulations with 20 CP units, CP1 still shows several free monomers together with some dimers, trimers, and a longer SCPN; CP2 forms 3 hexamers, together with some dimers and trimers; and CP3 forms mostly trimers, together with a pentamer and several dimers (Fig. S9). The simulations with 80 peptide units (Fig. S12) are more similar to those with 40 CPs. CP1 and CP3 form relatively small aggregates (mostly dimers and trimers) while CP2 forms predominantly tetramers, pentamers, hexamers and even heptamers. In the simulations containing 160 peptides, a large number of CP1 and CP3 remain as monomers at the end of the corresponding trajectories while CP2 tends to form SCPNs of 4 to 8 units (Fig. S14). For the simulations with the original MARTINI model using 40 CP units, CP2 and CP3 form medium sized SCPNs, while the behaviour of CP1 is more similar to that observed with the new parameterization (Fig. S11). The main assemblies observed in the trajectories with 20 CPs are again dimers and trimers for CP1, CP2 leads to one trimer, 2 pentamers and even one heptamer, and CP3 forms an octamer together with trimers and dimers (Fig. S10). When using 80 CP units, dimers and trimers are mostly formed for CP1 and CP3, while the assemblies with CP2 are larger, with SCPNs of 6, 9 and even 11 molecular units, and with no dimers and trimers at all (Fig. S13). Finally, the simulations with 160 units of CP1 and CP3 are again similar to those obtained with their respective studies in the new parametrization, being monomers, dimers and trimers the dominant species. CP2 changes its tendency a bit as the nanotubes formed are not as large as with the other concentrations, now 4-subunit SCPNs are predominant followed by dimers and trimers (Fig. S15). 3.2.1.2. Simulations in the presence of DMPC:DMPG 3:1 bilayers (M2). The behaviour of the peptides in this model of bacteria membrane is quantitatively different from that in the simulations with pure DMPC lipid bilayer (M1) but the general qualitative features persist. Again, for the simulations using 40 CP units with the MA(R/S)TINI parameterization, a significant number of CP1 and CP3 monomers survive as independent molecules while all the CP2 molecules aggregate forming medium size SCPNs (Fig. 4). Fig. 4. Number of CPs in clusters of different size (indicated by the colours in the legend) as a function of time, for simulations with 40 units of CP1,CP2 and CP3, in presence of the three membrane models, using the new MA(R/S)TINI CG parameterization developed in this work. Note that the initial value (at Time = 0) corresponding to the black line (representing the number of monomers) is always 40 for all the plots, while the rest of the lines start at zero. The exchange of CPs between clusters of different size can be clearly seen in the behaviour of the different curves. The results obtained for the rest of the free MD simulations with different number of CP units are in Figs S9-S15. A. Cabezón, M. Calvelo, J.R. Granja et al. Journal of Colloid and Interface Science 642 (2023) 84–99 91 The main quantitative difference is apparent for CP1, which shows a significantly larger number of monomers in the presence of this membrane than in the presence of M1: almost half of the peptides remain stable as independent structures and the rest of molecules form mostly dimers. The differences in the behaviour of CP2 and CP3 between the two membrane models are less clear. The trends of the simulations with 20 CP molecules are similar (Fig. S9), or even more marked than those using 40 peptide units. More than half of CP1 remain as monomers after 20 l s while no isolated CP2 molecules and just one dimer was found in the corresponding simulation and an intermediate behaviour is observed for CP3, with 6–7 stable monomers and a similar number of peptides assembled in dimers, trimers and one tetramer. Again, the dominant species are CP1 monomers when using 80 peptide units, medium size SCPNs for CP2 and a significant presence of monomers and mostly dimers for CP3 (Fig. S12). Monomers, dimers and trimers dominate the simulations of 160 CP1 and CP3 while SCPNs with 3 to 5 units are more likely for CP2 (Fig. S14). The simulations using the original MARTINI force field seem to favour the aggregation kinetics, compared to the results obtained using MA(R/S)TINI. The final number of CP1 monomers after 20 l s in the simulation with 40 peptide units is lower than that with the new parameterization, although most of the aggregates are small dimers and trimers (Fig. S11). A similar behaviour is observed for CP3, although for this peptide the SCPNs are a bit longer. CP2 exhibits a faster kinetics, reaching stable tetramers and pentamers just after 3 l s. Reducing the number of components to 20 in the simulation boxes, the three systems conserve the trends: CP1 forms small aggregates and less monomers than in the simulations with the new parameterization; CP2 forms longer aggregates of up to 9 subunits and, surprisingly, CP3 leads to an aggregate of 16 peptides that includes most of the CPs in the simulation box (Fig. S10). Differences between both parameterizations are smaller when the concentration of peptides increases up to 80 peptide molecules (Fig. S13). Slightly less monomers and more CP1 dimers are formed with MARTINI than with MA(R/S)TINI, 2 to 5 units length SCPNs are quickly formed with CP2, and CP3 is mainly distributed in monomers and dimers, with 3 trimers, 6 tetramers and even one hexamer at the end of the trajectory. For the simulations with 160 peptides (Fig. S15), the final distribution of structures is similar for CP1 and CP3, with most of the peptides as monomers, dimers and trimers, while CP2 aggregates much faster and forms significantly larger aggregates of up to 9 molecular units. 3.2.1.3. Simulations in the presence of DMPE:DMPG 1:9 bilayers (M3). The simulations with this lipid bilayer are characterized by a larger presence of CP1 monomers for both parameterizations of the force field, compared to the simulations with the other two membrane models. This is likely due to the electrostatic attraction between this cationic peptide and the DMPE:DMPG 1:9 bilayer (M3), which contains a high density of anionic lipids. Thus, the CP1 structures directly join the membrane as monomers, making difficult their aggregation. Typically, half of the CP1 structures remain as monomers throughout all the simulations, regardless the concentration and even of the force field parameterization, while the rest of the peptide units form mostly dimers and some trimers, also interacting with the membrane (Fig. 4 and S9-S15). CP2 keeps its ability to aggregate very quickly compared to CP1 and CP3, forming relatively large SCPNs of up to 8–10 peptide units for both parameterizations and for all the concentrations. Again, a key feature of CP3 is the persistency of monomers even at large concentrations of the peptide for both parameterizations. This behaviour is a constant for both CP1 and CP3, although to a lower extent for CP3. Differences in the kinetic behaviour between both force field parameterizations are not evident for any of the three CP sequences. The general trends are also similar, although quantitative differences can be observed. The behaviour of CP1 in the simulation with 20 peptide units is practically the same for both parameterizations (Fig. S9-S10) while for the rest of concentrations the ratio of CPs forming dimers to the number of monomers is higher for the original MARTINI parameterization than for MA(R/ S)TINI (Fig. 4 and S11-S15). CP3 also exhibits a persistent higher number of independent peptide units in the simulations with the new parameterization of the force field for all the studied concentrations. CP2, in agreement with the simulations in the other two membrane models, does not show free monomers at the end of the trajectories but the final aggregates are consistently larger at all the concentrations for the simulations with the original MARTINI force field. Overall, the kinetic analysis reveals interesting aspects on the competition between the adsorption of the peptides by the lipid bilayers and their self-assembly process. It is shown that the selfassembly kinetics is significantly faster for CP2 in practically all cases. This is probably connected to a specific interaction mechanism exhibited by this peptide (see next section). The selfassembly of CPs and the interaction with the lipid bilayers are seriously affected by the electrostatic attraction between different molecular entities (CP monomers, CP aggregates and membrane models). This is very clear for the most cationic peptide (CP1), which exhibits a progressively lower tendency to self-aggregate when increasing the negative charge density of the membrane, thus favouring its adsorption. It can be shown that the hydrophobic residues (W and L) of all the CPs in the SCPN are oriented towards the same side of the CP assemblies using both MA(R/S) TINI and MARTINI (Fig. S16). When a lipid membrane is present, these hydrophobic residues interact with the lipid tails of the bilayer, while the ionic residues (R and K) are exposed to the polar solvent. One difference between these two forcefields is that with the new parametrization, the formation of an amphipathic SCPN, i.e., with polar and nonpolar residues aligned with themselves, seems to be less likely in the absence of the membrane. This means that the relative orientation of the residues is more random with MA(R/S)TINI than with MARTINI in aqueous media. Additionally, the length of the final SCPNs is also consistently larger for the neutral CP2 than for the charged CP1 and, to a lower extent for the polar CP3. Globally, no significant quantitative differences can be observed between MARTINI and MA(R/S)TINI in terms of the number of CP clusters formed during the simulated time (Fig. 5). 3.2.2. Final self-assembly patterns The kinetic quantitative results are clearly connected to the final self-assembly patterns achieved by each system (Fig. 6). When using the chirality aware MA(R/S)TINI force field in presence of the mammalian membrane model (M1), the cationic CP1 and the non-ionic but polar CP3 tend to aggregate in the aqueous solution more than to interact with the neutral lipid bilayer, while CP2 clearly exhibits two different interaction mechanisms: formation of SCPNs throughout the lipid membranes and formation of twodimensional lattices in water solution. This behaviour persists for all the peptide concentrations. In contrast, CP1 is strongly attracted for both bacteria membrane models (M2 and M3) favoured by electrostatic interactions. When the concentration increases, CP1 seriously perturbs the structure of such membrane models, mainly that of M3, indicating the specificity of this peptide towards membranes of different composition and, in particular, charge density. These results suggest the interaction mechanism for the toxicity of CP1 towards bacteria. The simulations using the original MARTINI force field provide similar results to those obtained with the MA(R/S)TINI parameterization (Fig. S17). At low concentrations A. Cabezón, M. Calvelo, J.R. Granja et al. Journal of Colloid and Interface Science 642 (2023) 84–99 92 [28] C.A. López, A.J. Rzepiela, A.H. de Vries, L. Dijkhuizen, P.H. Hünenberger, S.J. Marrink, Martini Coarse-Grained Force Field: Extension to Carbohydrates, J. Chem. Theory Comput. 5 (12) (2009) 3195–3210, https://doi.org/10.1021/ ct900313w. [29] L. Monticelli, S.K. Kandasamy, X. Periole, R.G. Larson, D.P. Tieleman, S.J. Marrink, The MARTINI Coarse-Grained Force Field: Extension to Proteins, J. Chem. Theory Comput. 4 (5) (2008) 819–834, https://doi.org/10.1021/ ct700324x. [30] P.C.T. Souza, R. Alessandri, J. Barnoud, S. Thallmair, I. Faustino, F. Grünewald, I. Patmanidis, H. Abdizadeh, B.M.H. Bruininks, T.A. Wassenaar, P.C. Kroon, J. Melcr, V. Nieto, V. Corradi, H.M. Khan, J. Doman ´ski, M. Javanainen, H. MartinezSeara, N. Reuter, R.B. Best, I. Vattulainen, L. Monticelli, X. Periole, D.P. Tieleman, A.H. de Vries, S.J. Marrink, Martini 3: A General Purpose Force Field for CoarseGrained Molecular Dynamics, Nat. Methods 18 (4) (2021) 382–388, https:// doi.org/10.1038/s41592-021-01098-3. [31] D.H. De Jong, G. Singh, W.F.D. Bennett, C. Arnarez, T.A. Wassenaar, L.V. Schäfer, X. Periole, D.P. Tieleman, S.J. Marrink, Improved Parameters for the Martini Coarse-Grained Protein Force Field, J. Chem. Theory Comput. 9 (1) (2013) 687– 697, https://doi.org/10.1021/ct300646g. [32] S.J. Marrink, L. Monticelli, M.N. Melo, R. Alessandri, D.P. Tieleman, P.C.T. Souza, Two Decades of Martini: Better Beads, Broader Scope, Wiley Interdiscip. Rev.: Comput. Mol. Sci. John Wiley & Sons, Ltd (2022) e1620. [33] S.J. Marrink, H.J. Risselada, S. Yefimov, D.P. Tieleman, A.H. De Vries, The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations, J. Phys. Chem. B 111 (27) (2007) 7812–7824, https://doi.org/10.1021/jp071097f. [34] C. Wang, L. Guan, D. Danovich, S. Shaik, Y. Mo, The Origins of the Directionality of Noncovalent Intermolecular Interactions, J. Comput. Chem. 37 (1) (2016) 34–45, https://doi.org/10.1002/jcc.23946. [35] P. Tiangtrong, N. Thamwattana, D. Baowan, Modelling Water Molecules inside Cyclic Peptide Nanotubes, Appl. Nanosci. 6 (3) (2016) 345–357, https://doi. org/10.1007/s13204-015-0436-4. [36] M. Tarek, B. Maigret, C. Chipot, Molecular Dynamics Investigation of an Oriented Cyclic Peptide Nanotube in DMPC Bilayers, Biophys. J. 85 (4) (2003) 2287–2298, https://doi.org/10.1016/S0006-3495(03)74653-0. [37] A. Khalfa, M. Tarek, On the Antibacterial Action of Cyclic Peptides: Insights from Coarse-Grained MD Simulations, J. Phys. Chem. B 114 (8) (2010) 2676– 2684, https://doi.org/10.1021/jp9064196. [38] S. Fernandez-Lopez, H.S. Kim, E.C. Choi, M. Delgado, J.R. Granja, A. Khasanov, K. Kraehenbuehl, G. Long, D.A. Weinberger, K.M. Wilcoxen, M.R. Ghadiri, Antibacterial Agents Based on the Cyclic D,La -Peptide Architecture, Nature 412 (6845) (2001) 452–455, https://doi.org/10.1038/35086601. [39] G. Van Meer, D.R. Voelker, G.W. Feigenson, Membrane Lipids: Where They Are and How They Behave, Nat. Rev. Mol. Cell Bio. (2008) 112–124, https://doi.org/ 10.1038/nrm2330. [40] G. Shahane, W. Ding, M. Palaiokostas, M. Orsi, Physical Properties of Model Biological Lipid Bilayers: Insights from All-Atom Molecular Dynamics Simulations, J. Mol. Model. 25 (3) (2019) 1–13, https://doi.org/10.1007/ s00894-019-3964-0. [41] K. Murzyn, T. Róg, M. Pasenkiewicz-Gierula, PhosphatidylethanolaminePhosphatidylglycerol Bilayer as a Model of the Inner Bacterial Membrane, Biophys. J. 88 (2) (2005) 1091–1103, https://doi.org/10.1529/ biophysj.104.048835. [42] M.Z. Tien, D.K. Sydykova, A.G. Meyer, C.O. Wilke, Peptidebuilder: A Simple Python Library to Generate Model Peptides, PeerJ 2013 (1) (2013) e80. [43] P. Gkeka, L. Sarkisov, Interactions of Phospholipid Bilayers with Several Classes of Amphiphilic a -Helical Peptides: Insights from Coarse-Grained Molecular Dynamics Simulations, J. Phys. Chem. B 114 (2010) 826–839, https://doi.org/ 10.1021/jp908320b. [44] M.J. Abraham, T. Murtola, R. Schulz, S. Páll, J.C. Smith, B. Hess, E. Lindah, Gromacs: High Performance Molecular Simulations through Multi-Level Parallelism from Laptops to Supercomputers, SoftwareX 1–2 (2015) 19–25, https://doi.org/10.1016/j.softx.2015.06.001. [45] Y. Qi, H.I. Ingólfsson, X. Cheng, J. Lee, S.J. Marrink, W. Im, CHARMM-GUI Martini Maker for Coarse-Grained Simulations with the Martini Force Field, J. Chem. Theory Comput. 11 (9) (2015) 4486–4494, https://doi.org/10.1021/acs. jctc.5b00513. [46] H.J.C. Berendsen, J.P.M. van Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, Molecular Dynamics with Coupling to an External Bath, J. Chem. Phys. 81 (8) (1984) 3684–3690. [47] M. Parrinello, A. Rahman, Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method, J. Appl. Phys. 52 (12) (1981) 7182–7190, https:// doi.org/10.1063/1.328693. [48] G. Bussi, D. Donadio, M. Parrinello, Canonical Sampling through Velocity Rescaling, J. Chem. Phys. 126 (1) (2007), https://doi.org/10.1063/1.2408420 014101. [49] B. Hess, H. Bekker, H.J.C. Berendsen, J.G.E.M. Fraaije, LINCS: A Linear Constraint Solver for Molecular Simulations, J. Comput. Chem. 18 (18) (1997) 1463–1472, https://doi.org/10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0. CO;2-H. [50] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, A Smooth Particle Mesh Ewald Method, J. Chem. Phys. 103 (19) (1995) 8577– 8593, https://doi.org/10.1063/1.470117. [51] S.O. Yesylevskyy, L.V. Schäfer, D. Sengupta, S.J. Marrink, Polarizable Water Model for the Coarse-Grained MARTINI Force Field, PLoS Comput. Biol. 6 (6) (2010) 1–17, https://doi.org/10.1371/journal.pcbi.1000810. [52] A. Lindahl, S.V. de Hess, GROMACS 2021 Manual. Zenodo. GROMACS 2021 Manual. Zenodo 2021. [53] M. Bonomi, G. Bussi, C. Camilloni, G.A. Tribello, P. Banáš, A. Barducci, M. Bernetti, P.G. Bolhuis, S. Bottaro, D. Branduardi, R. Capelli, P. Carloni, M. Ceriotti, A. Cesari, H. Chen, W. Chen, F. Colizzi, S. De, M. De La Pierre, D. Donadio, V. Drobot, B. Ensing, A.L. Ferguson, M. Filizola, J.S. Fraser, H. Fu, P. Gasparotto, F.L. Gervasio, F. Giberti, A. Gil-Ley, T. Giorgino, G.T. Heller, G.M. Hocky, M. Iannuzzi, M. Invernizzi, K.E. Jelfs, A. Jussupow, E. Kirilin, A. Laio, V. Limongelli, K. Lindorff-Larsen, T. Löhr, F. Marinelli, L. Martin-Samos, M. Masetti, R. Meyer, A. Michaelides, C. Molteni, T. Morishita, M. Nava, C. Paissoni, E. Papaleo, M. Parrinello, J. Pfaendtner, P. Piaggi, G.M. Piccini, A. Pietropaolo, F. Pietrucci, S. Pipolo, D. Provasi, D. Quigley, P. Raiteri, S. Raniolo, J. Rydzewski, M. Salvalaglio, G.C. Sosso, V. Spiwok, J. Šponer, D.W.H. Swenson, P. Tiwary, O. Valsson, M. Vendruscolo, G.A. Voth, A. White, Promoting Transparency and Reproducibility in Enhanced Molecular Simulations, Nat. Methods (2019) 670– 673, https://doi.org/10.1038/s41592-019-0506-8. [54] W. Humphrey, A. Dalke, K. Schulten, VMD: Visual Molecular Dynamics, J. Mol. Graph. 14 (1) (1996) 33–38, https://doi.org/10.1016/0263-7855(96)00018-5. [55] N. Michaud-Agrawal, E.J. Denning, T.B. Woolf, O. Beckstein, MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations, J. Comput. Chem. 32 (10) (2011) 2319–2327, https://doi.org/10.1002/jcc.21787. [56] C.R. Harris, K.J. Millman, S.J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N.J. Smith, R. Kern, M. Picus, S. Hoyer, M.H. van Kerkwijk, M. Brett, A. Haldane, J.F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, T.E. Oliphant, Array Programming with NumPy, Nature (2020) 357–362, https://doi.org/10.1038/s41586-020-2649-2. [57] J.D. Hunter, Matplotlib: A 2D Graphics Environment, Comput. Sci. Eng. 9 (3) (2007) 90–95, https://doi.org/10.1109/MCSE.2007.55. [58] This type of peptides is very poorly soluble in water and tends to precipitate as soon as CP solutions in DMSO are added to aqueous solutions containing liposomes, so their behavior as transmembrane channel formers is unknown (unpublished results). [59] I. Insua, J. Montenegro, 1D to 2D Self Assembly of Cyclic Peptides, J. Am. Chem. Soc. 142 (1) (2020) 300–307, https://doi.org/10.1021/jacs.9b10582. A. Cabezón, M. Calvelo, J.R. Granja et al. Journal of Colloid and Interface Science 642 (2023) 84–99 99