Literature DB >> 33864790

Enhanced translocation of amphiphilic peptides across membranes by transmembrane proteins.

Ladislav Bartoš1, Ivo Kabelka1, Robert Vácha2.   

Abstract

Cell membranes are phospholipid bilayers with a large number of embedded transmembrane proteins. Some of these proteins, such as scramblases, have properties that facilitate lipid flip-flop from one membrane leaflet to another. Scramblases and similar transmembrane proteins could also affect the translocation of other amphiphilic molecules, including cell-penetrating or antimicrobial peptides. We studied the effect of transmembrane proteins on the translocation of amphiphilic peptides through the membrane. Using two very different models, we consistently demonstrate that transmembrane proteins with a hydrophilic patch enhance the translocation of amphiphilic peptides by stabilizing the peptide in the membrane. Moreover, there is an optimum amphiphilicity because the peptide could become overstabilized in the transmembrane state, in which the peptide-protein dissociation is hampered, limiting the peptide translocation. The presence of scramblases and other proteins with similar properties could be exploited for more efficient transport into cells. The described principles could also be utilized in the design of a drug-delivery system by the addition of a translocation-enhancing peptide that would integrate into the membrane.
Copyright © 2021 Biophysical Society. Published by Elsevier Inc. All rights reserved.

Entities:  

Mesh:

Substances:

Year:  2021        PMID: 33864790      PMCID: PMC8390799          DOI: 10.1016/j.bpj.2021.04.005

Source DB:  PubMed          Journal:  Biophys J        ISSN: 0006-3495            Impact factor:   3.699


Significance

Cells are separated from the external environment by a selectively permeable cytoplasmic membrane. Peptides with the right properties can spontaneously cross this protective barrier and thus act as drug carriers or therapeutics themselves. However, the right properties remain enigmatic. Here, we show that amphiphilic peptides without such properties could also spontaneously enter the cell with the passive help of transmembrane proteins present in the membrane. Translocation-enhancing membrane proteins create a passage that can be utilized for the facilitated permeation of amphiphilic molecules into cells. Moreover, peptides with properties similar to the studied transmembrane proteins could be used in synergistic mixtures to enhance the translocation.

Introduction

Cell membranes regulate the exchange of matter between cells and their environment. The membranes are semipermeable, enabling some molecules to spontaneously pass through this protective barrier. These molecules are typically small and uncharged. However, larger amphiphilic molecules, including peptides, were also reported to spontaneously translocate through membranes (1, 2, 3). Such peptides include antimicrobial peptides (AMPs) and cell-penetrating peptides (CPPs), which are promising drugs and drug carriers, respectively (4, 5, 6, 7). Despite recent progress, it remains unclear what the right physicochemical properties are to enable spontaneous peptide translocation into the cells (1,3,8,9). Even peptides/molecules that are not able to spontaneously translocate across membranes can be internalized into a cell. One obvious way is a tightly regulated transporter-mediated uptake. However, a local decrease in membrane hydrophobic thickness can also facilitate a passive uptake by lowering the energetic barrier. Such a decrease can be caused by local membrane thinning due to the presence of shorter and oxidized lipids, membrane proteins, and/or other amphiphilic molecules (3,9). Integral membrane proteins with partly hydrophilic transmembrane domain(s), such as scramblases, are particularly interesting. In such proteins, the presence of hydrophilic residues in a membrane locally increases polarity, which can enhance lipid flip-flop, a translocation of lipids between membrane leaflets (10, 11, 12, 13). We hypothesize that scramblases and similar proteins can also enhance the translocation of other amphiphilic molecules, including peptides that would not otherwise be able to spontaneously translocate across membranes (see Fig. 1). A similar mechanism could be utilized in peptide mixtures in which one peptide would embed itself into the membrane and facilitate the translocation of other peptides. Therefore, our results could be used in the rational design of novel AMPs/CPPs and their mixtures.
Figure 1

Hypothesis of peptide translocation enhanced by transmembrane protein/peptide demonstrated on the free energy profiles of peptide translocation as functions of a distance between the peptide and membrane centers of mass. The translocating peptide (TLP) has a high barrier for passing through the membrane individually (red). When a transmembrane protein/peptide (MP) is present in the membrane, the translocation of TLP is enhanced due to membrane disruption caused by the MP and due to the stabilizing interaction between the TLP and MP (green). To see this figure in color, go online.

Hypothesis of peptide translocation enhanced by transmembrane protein/peptide demonstrated on the free energy profiles of peptide translocation as functions of a distance between the peptide and membrane centers of mass. The translocating peptide (TLP) has a high barrier for passing through the membrane individually (red). When a transmembrane protein/peptide (MP) is present in the membrane, the translocation of TLP is enhanced due to membrane disruption caused by the MP and due to the stabilizing interaction between the TLP and MP (green). To see this figure in color, go online. Here, we employed two different models to calculate the free energy profiles of peptide translocation across phospholipid membranes in the presence of a transmembrane protein/peptide. We focused on α-helical peptides with various hydrophilicities and analyzed their effect on the free energy profiles of translocation with a focus on the process limiting barriers. By comparing our phenomenological model with a Martini coarse-grained model, we concluded that partly hydrophilic transmembrane proteins/peptides can act as translocation enhancers that facilitate the translocation of amphiphilic peptides into the cells.

Materials and methods

Phenomenological model

The simulations were performed with an in-house-developed, freely available software, SC (https://github.com/robertvacha/SC), using the Metropolis Monte Carlo algorithm (14). For lipids, we employed an implicit-solvent Cooke-Deserno model (15), which describes each lipid as a chain of three spherical particles. The hydrophilic lipid headgroup is represented by a single, purely repulsive particle with a diameter of 0.95 nm. Both hydrophobic tails are fused together and modeled by two attractive particles with a diameter of 1.0 nm. This lipid model reproduces the elastic properties and phase transition of a lipid membrane (15). A phenomenological model of patchy spherocylinders (PSC) (16) was used to represent both translocating peptides (TLPs) and transmembrane proteins/peptides (MPs). This model was developed to describe amphiphilic helices and to be compatible with the lipid model by Cooke and Deserno (15). The entire peptide or protein is modeled by a single spherocylinder particle containing a hydrophobic (attractive) patch, which is defined by the angle of the cylinder wedge (see Fig. 2 A), and hydrophilic (purely repulsive) patch, which covers the rest of the surface. The hydrophobic patch can either span the entire length of the spherocylinder, including the hemispherical endcaps, or it can only span the cylindrical segment of the particle. In the latter case, the hemispherical endcaps are purely repulsive, representing hydrophilic termini of the peptide/protein. The first type of spherocylinder is denoted PSC with attractive endcaps (PSC-AE), and the second type is called PSC with nonattractive endcaps (PSC-NE) (see Fig. 2 B).
Figure 2

Illustration of the employed models and sequences. In the phenomenological model, an α-helical amphiphilic peptide is represented by a patchy spherocylinder (PSC) with a certain size of the hydrophobic patch (gray) given by its angular wedge (A) and with a specific character of the endcaps (B): hydrophilic, NE, or partially hydrophobic, AE. (C) and (D) show helical wheel representations of sequences of translocating peptides (TLPs) and transmembrane proteins/peptides (MPs), respectively, simulated with Martini force field. Peptide/protein names are based on the SAG terminal sequence (only for MPs), polyleucine core, and the number of serine residues forming the hydrophilic part of the peptide. To see this figure in color, go online.

Illustration of the employed models and sequences. In the phenomenological model, an α-helical amphiphilic peptide is represented by a patchy spherocylinder (PSC) with a certain size of the hydrophobic patch (gray) given by its angular wedge (A) and with a specific character of the endcaps (B): hydrophilic, NE, or partially hydrophobic, AE. (C) and (D) show helical wheel representations of sequences of translocating peptides (TLPs) and transmembrane proteins/peptides (MPs), respectively, simulated with Martini force field. Peptide/protein names are based on the SAG terminal sequence (only for MPs), polyleucine core, and the number of serine residues forming the hydrophilic part of the peptide. To see this figure in color, go online. Each system was composed of 500 lipid molecules assembled in a bilayer in a rectangular box with dimensions of roughly 17×17×35 nm. The membrane normal was oriented along the z axis. The first spherocylinder, representing the MP, was embedded into the center of the membrane in an orientation perpendicular to the membrane plane. The second spherocylinder, representing the TLP, was placed slightly above the membrane surface in an orientation parallel to the membrane plane and close to the MP. All simulations were performed in an NPT ensemble with a reduced temperature of 1.0 kT/ε and reduced pressure in the xy plane equal to 0.0 pσ 3/ε, which corresponds to zero tension in the membrane. σ is a unit of length that roughly corresponds to 1 nm (15). Each spherocylinder had a diameter of 1 nm, roughly matching the size of an ideal α-helix. The length of helices/spherocylinders was 4 nm, corresponding to the width of the hydrophobic core of the membrane. All transmembrane helices (MPs) were of the NE type, i.e., with hydrophilic ends. For TLPs, we studied both PSC-AE and PSC-NE. Hydrophobic patch sizes of 180, 270, and 360° were used for MPs; TLPs had hydrophobic patch sizes of 90, 135, 180, and 270°. Every combination of MP and TLP was simulated. We also prepared one reference system without any MP for each TLP. The MP was fixed in the xy position, and we restrained its movement on the z axis to the membrane center of mass using a flat-bottom potential:with the force constant k being equal to 10ε, d corresponding to the distance between the MP center of mass and membrane center of mass, and dref = 0.5 nm. To enhance the sampling, the displacement of the TLP in the xy plane was restrained to the neighborhood of the MP with a flat-bottom potential. In this case, d corresponded to the xy distance between the TLP center of mass and MP center of mass, and the reference value dref was 4 nm. The TLP could translocate either directly along MP or individually through the membrane but could not diffuse far away.

Free energy calculations

The Wang-Landau method (17) was employed for free energy calculations. To describe peptide translocation, we used a previously designed two-dimensional collective variable (CV) (18). The first CV dimension was the oriented distance of the TLP center of mass from the membrane center of mass in the z axis. The second CV dimension was the orientation of the TLP in z axis, defined as the cosine of the angle between the main peptide axis and the z axis, . The range of the first CV dimension was −5.7 to 5.7 nm with a bin width of 0.1 nm, and the range of the second CV dimension was 0–1 with a bin width of 0.0125. The initial modification factor f for the Wang-Landau method was set to exp(10-6) and was decreased to whenever the following two convergence conditions were met: 1) each histogram bin had at least 1000 samples and 2) the maximal roughness of the histogram defined as kT × log(H/H) was lower than 10−4 kT, where k is the Boltzmann constant, T is the thermodynamic temperature, and H is the number of samples in the specified histogram bin. Once the modification factor f decreased, the histogram was emptied. The f factor stopped decreasing once it became lower than exp(10-7). Subsequently, we carried out simulations without modifying the f factor (with the detailed balance) and collected histograms in CVs. These histograms were used to ensure the quality of sampling and to refine the free energy surfaces. In one-dimensional free energy profiles, the orientations of the TLP were averaged out. Free energy was calculated in kT and then converted to kJ/mol using T = 310 K.

Calculations of lipid defect

To quantify the effect of MP on membrane structure, we simulated one system for each of the spherocylinders 180, 210, 240, 270, and 360° PSC-NE embedded in the membrane. Each system was prepared in the same way as described above, but no TLP was placed in the system, and the Wang-Landau method was not used. We then simulated each system for 20 million Monte Carlo sweeps, taking a snapshot every 5000 sweeps. In each snapshot, we calculated the density of lipid heads depending on their z distance from the membrane center of mass. The “lipid defect” was then calculated as the average number of lipid heads located both within 2.5 nm from the MP in the xy plane and within 1.5 nm from the membrane center of mass on the z axis.

Martini simulations

All our molecular dynamics simulations were performed with Gromacs package version 5.1.4 (19). We employed coarse-grained Martini force field version 2.2 (20, 21, 22). A recently proposed correction, which scaled down the Lennard-Jones interactions between peptide and protein beads by 10%, was applied to prevent unrealistically strong protein-protein interactions (23). Peptides and proteins were prepared using MODELER version 9.11 (24) and then coarse-grained using the martinize.py script version 2.6 (https://github.com/cgmartini/martinize.py). Each peptide/protein was then minimized in vacuum using the steepest-descent algorithm until the maximal force on any atom was lower than 100 kJ mol−1 nm−1. Hydrophobicity of each peptide/protein was calculated using HeliQuest ComputParams.py script version 3 (https://heliquest.ipmc.cnrs.fr/cgi-bin/ComputParams.py). The employed hydrophobicity scale is based on the partitioning of N-acetyl-amino-acid amides in octanol (25). Each system consisted of a bacterial-membrane mimic composed of 96 molecules of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine and 32 molecules of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoglycerol in each leaflet (256 lipids in the whole membrane). The bilayer was prepared using the CHARMM-GUI web interface (26) in the xy plane with its normal oriented along the z axis. For insertion simulations, the system had an approximate size of 9×9×11 nm and was solvated with roughly 4400 water beads. For adsorption simulations, the system was larger, at 9×9×22 nm, and it was solvated with roughly 11,700 water beads. The 23-amino-acid-long α-helical transmembrane protein/peptide (MP) was placed into the center of the membrane in an orientation perpendicular to the membrane plane. The 21-amino-acid-long α-helical TLP was initially positioned on the surface of the membrane in an orientation parallel to the membrane plane and with its hydrophobic patch oriented toward the membrane. For adsorption simulations, another identical TLP was placed on the surface of the opposite membrane leaflet to avoid any discrepancies caused by changes in the leaflet tension. For each translocating peptide, a reference system without MP was also simulated. The names and sequences of all simulated TLPs and MPs are shown in Tables 1 and 2, respectively, and graphically depicted in Fig. 2, C and D.
Table 1

Names, sequences, and hydrophobicities of simulated TLPs

NameSequenceHydrophobicity
LS5LSLLLLLLSLLLSLLSLLLSL-NH21.286
LS9LSSLLSLLSSLLSLLSSLLSL-NH20.954
LS13SSSLSSLLSSLSSLLSSLSSL-NH20.623

Names are derived from polyleucine sequence with the specific number of serine residues.

Table 2

Names, sequences, and hydrophobicities of simulated MPs

NameSequenceHydrophobicity
SAGLSAGLLLLLLLLLLLLLLLLLGAS1.280
SAGLS5SAGLSLLLLLLSLLLLLLSLGAS1.053
SAGLS7SAGLSLLSLLLSLLLSLLSLGAS0.902
SAGLS9SAGLSLLSSLLSLLSSLLSLGAS0.750
SAGLS11SAGSSLLSSLLSLLSSLLSSGAS0.599
SAGLS13SAGSSLLSSLSSSLSSLLSSGAS0.448
SAGLS15SAGSSSLSSLSSSLSSLSSSGAS0.297

Names are derived from SAG termini and polyleucine sequence with the specific number of serine residues.

Names, sequences, and hydrophobicities of simulated TLPs Names are derived from polyleucine sequence with the specific number of serine residues. Names, sequences, and hydrophobicities of simulated MPs Names are derived from SAG termini and polyleucine sequence with the specific number of serine residues. Finally, Na+ and Cl− ions were added to the system at physiological concentration of 154 mM with an excess of ions to neutralize the system. Subsequently, each system was minimized using the steepest-descent algorithm with maximal force tolerance of 100 kJ mol−1 nm−1. Equilibration was performed in five stages with increasing simulation time steps and various overall lengths: 1) dt = 2 fs and t = 0.5 ns, 2) dt = 5 fs and t = 1.25 ns, 3) dt = 10 fs and t = 1 ns, 4) dt = 20 fs and t =30 ns, and 5) dt = 20 fs and t = 15 ns. In all but the last stage, position restraints in the form of a harmonic potential with a force constant of 1000 kJ mol−1 nm−2 were applied to all backbone bead coordinates of both MP and TLP(s). In stage 5 and all subsequent simulations, TLP was forced to stay close to the MP using a flat-bottom potential with a force constant of 500 kJ mol−1 nm−2. The reference distance between the center of mass of the TLP backbone and the center of mass of the MP backbone was set to 2.5 nm in the xy plane. In adsorption simulations, the TLPs on both membrane leaflets were restrained using a flat-bottom potential to remain in proximity to the MP. During the equilibration, another restraint was applied to the termini of both TLPs, which were held at a distance of 1.4 nm from the local membrane center of mass on the z axis using a harmonic potential with a force constant of 500 kJ mol−1 nm−2. This distance is slightly deeper than the equilibrium position of peptide termini, and the restraint was applied to ensure sufficient overlap between the free energy profiles of adjacent adsorption and insertion simulations. Equilibration was done in an NPT ensemble with the temperature being maintained at 310 K using a modified velocity-rescaling thermostat (27) with a coupling constant of 1 ps. Water with ions and membrane with the TLP and the MP were coupled to two separate thermal baths. Pressure was kept at 1 bar using Berendsen barostat (28) with coupling constant of 12 ps. Semi-isotropic pressure coupling was employed to independently scale the simulation box in the xy plane and on the z axis with a compressibility of 3×10−4 bar−1. A leap-frog algorithm was used to integrate the Newtonian equations of motion. Nonbonded interactions were cut off at 1.1 nm. Van der Waals potential was shifted to zero at cutoff distance. The relative dielectric constant was set to 15. After the equilibration, a short 100-ns molecular dynamics was carried out. For this run, as well as for the subsequent pulling and umbrella sampling simulations, a Berendsen barostat was replaced with a Parinello-Rahman barostat (29,30). All other simulation settings were the same as in stage 5 of the equilibration. To enhance the sampling of the configuration space, we employed an umbrella sampling method (31,32). As previously (8), we divided the process of peptide translocation into four separate processes: C-terminus adsorption, C-terminus insertion, N-terminus insertion, and N-terminus adsorption. Because the membrane and the MP were both symmetric, the free energy profile of the full translocation pathway could be obtained by joining the free energy profiles of individual insertions and adsorptions. The CV, chosen to describe each of the insertions/adsorptions, was the oriented distance between the terminus center of mass and the local membrane center of mass on the z axis. The terminus was defined as the first (N) or the last (C) three backbone beads of the translocating peptide. The local membrane center of mass was calculated from the positions of lipid beads localized in a cylinder with the radius of 2.0 nm and its principal axis going through the TLP terminus. For insertion simulations, the initial configurations for umbrella sampling were generated by pulling a single TLP terminus across the membrane. The terminus was pulled for approximately 1 μs with a pulling rate of 4.2 nm μs−1 and an initial reference distance of 2.1 nm using a harmonic potential with a force constant of 8000 kJ mol−1 nm−2. For adsorption simulations, the configurations for umbrella sampling were obtained by pulling the corresponding termini of two TLPs positioned on opposite leaflets in opposite directions. Termini were pulled for ∼400 ns with a pulling rate of 15 nm μs−1 and the initial reference distance of ±1.4 nm using a harmonic potential with a force constant of 8000 kJ mol−1 nm−2. The insertion pulling trajectory was then split into 64 nonuniformly distributed sampling windows (see Table S1), and the adsorption pulling trajectory was split into 56 sampling windows, as shown in Table S2. After a short 30-ns equilibration, each insertion window was sampled for at least 1 μs, and each adsorption window was sampled for 500 ns. The umbrella force constant was set to 8000 kJ mol−1 nm−2 for windows near the bilayer center and to 3000 kJ mol−1 nm−2 for windows near the membrane surface. For some systems, a force constant of 5000 kJ mol−1 nm−2 was applied in windows bordering these two regions. For adsorption windows, we used force constants between 1000 and 1800 kJ mol−1 nm−2. The full free energy profile for each N-/C-terminus insertion or N-/C-terminus adsorption was obtained from a set of umbrella windows using the weighted histogram analysis method (33,34) as implemented in g_wham (35). For the adsorption simulations, each set of umbrella sampling simulations led to two free energy profiles, one for each peptide. The final profile for N-/C-terminus adsorption was obtained by taking an average of them. The free energy profile was checked for artifacts from membrane size by running an additional simulation with the membrane having double the number of lipids (512 lipids) (see Fig. S1).

Calculations of water defect

To characterize the effect of the presence of the MP on membrane structure, we simulated one system for each MP in which single MP was placed into the membrane but no TLP was present. The system was prepared, minimized, and equilibrated in the same way as described above, but the molecular dynamics stage was extended to 1 μs. We analyzed the molecular dynamics trajectory, taking a snapshot every 200 ps and calculating the density of water beads depending on their z distance from the membrane center of mass. The water defect around the MP was calculated as the average number of water beads located within 2.5 nm from the MP in the xy plane and within 2 nm from the membrane center of mass on the z axis. To further characterize the role of the water defect in peptide translocation, we calculated the free energy of insertion for the peptide LS9 in the presence of three different artificially created water defects that mimicked the presence of the MPs SAGL, SAGLS7, and SAGLS15. Water defects were created by restraining all beads of lipid tails using an inverted flat-bottom potential in the center of the membrane. The potential had a force constant of 8 (mimicking SAGL), 14.1 (mimicking SAGLS7), or 23 kJ mol−1 nm−2 (mimicking SAGLS15). The reference distance between the lipid beads and the center of the system in the xy plane was 1.0 nm. Similar to the translocation in the presence of the MP, the translocating peptide was restrained in its movement in the xy plane. Instead of the MP, the restraining potential was centered on the defect by using a flat-bottom potential with a reference distance of 2.5 nm and a force constant of 200 kJ mol−1 nm−2 applied to three backbone beads of residues in the center of the peptide sequence. The free energy of peptide insertion for these systems was calculated in the same way as described above.

Results

We calculated the free energy of translocation of an α-helical TLP across a phospholipid membrane containing an α-helical transmembrane protein/peptide (MP). The hydrophobicity of both the TLP and the MP was varied, but MP always had hydrophilic ends (NE); see the Materials and methods for details. The translocation process followed the same pathway in all our studied cases. First, the TLP adsorbed onto the surface of the membrane in an orientation parallel to the membrane plane. Then, the TLP inserted one of its termini into the hydrophobic core of the membrane while changing its orientation. In the transmembrane state, the TLP was oriented perpendicular to the membrane plane. In most systems, the TLP and MP were oriented with their hydrophilic patches toward each other during the translocation through the membrane. The free energy of the TLP was decreased the most in the transmembrane state, in which TLP and MP were aligned with the membrane normal (see Fig. S2 for simulation snapshots depicting the insertion mechanism). The exception was systems containing an MP with no hydrophilic patch (MP with 360° hydrophobic patch). In such systems, the free energy of the TLP slightly increased, and the TLP usually translocated through the membrane alone, not directly along the MP. Representative free energy profiles for translocation of a half-hydrophobic/half-hydrophilic TLP (180° PSC-NE) are shown in Fig. 3. The presence of the MP had a dramatic effect on the translocation. The TLP free energy in the membrane decreased with increasing hydrophilicity of the MP, i.e., the TLP was more stable in the membrane in the presence of a more hydrophilic MP. The same trend was also observed for other investigated TLPs, irrespective of the hydrophilicity of their termini (NE/AE). Highly hydrophobic TLPs were slightly less affected by the presence of the MP in the membrane than more hydrophilic translocating peptides. See Figs. S3–S10 for two-dimensional free energy surfaces of TLPs in various systems and Fig. S11 for one-dimensional free energy profiles with averaged-out peptide orientations.
Figure 3

Free energy profile of half-hydrophobic/half-hydrophilic peptide (180° PSC-NE) translocating in the presence of various transmembrane proteins/peptides (MPs) with different hydrophobicities (180, 270, and 360° PSC-NE). The free energy is shown as a function of the translocating peptide distance from the membrane center of mass. The black line corresponds to a reference system without any MPs. The translocation process is depicted in the illustration at the bottom. Hydrophobic patches of the peptides are shown in gray, and hydrophilic areas are colored green. The dark-cyan arrow shows the free energy barrier (ΔΔGBT) of the translocation profile for the TLP 180° PSC-NE in the presence of the MP 270° PSC-NE. ΔΔGBT is defined as the difference between the global maximum and minimum in the free energy profile. The error of profiles was estimated to be below 2.5 kJ/mol based on the difference between solution states (profile ends). To see this figure in color, go online.

Free energy profile of half-hydrophobic/half-hydrophilic peptide (180° PSC-NE) translocating in the presence of various transmembrane proteins/peptides (MPs) with different hydrophobicities (180, 270, and 360° PSC-NE). The free energy is shown as a function of the translocating peptide distance from the membrane center of mass. The black line corresponds to a reference system without any MPs. The translocation process is depicted in the illustration at the bottom. Hydrophobic patches of the peptides are shown in gray, and hydrophilic areas are colored green. The dark-cyan arrow shows the free energy barrier (ΔΔGBT) of the translocation profile for the TLP 180° PSC-NE in the presence of the MP 270° PSC-NE. ΔΔGBT is defined as the difference between the global maximum and minimum in the free energy profile. The error of profiles was estimated to be below 2.5 kJ/mol based on the difference between solution states (profile ends). To see this figure in color, go online. The efficiency of peptide translocation is mainly determined by the free energy barrier, ΔΔGBT, on the pathway. The height of the translocation barrier is the difference between the maximum and minimum, ΔGMAX−ΔGMIN, on the free energy profile. Free energy minima, maxima, and barriers, as well as the free energy of the transmembrane and (meta)stable adsorbed states, are summarized in Table 3 for all simulated systems. As can be seen, ΔΔGBT does not always decrease with increasing hydrophilicity of the embedded MP. This is caused by the fact that highly hydrophilic MPs can overstabilize the TLP in the membrane hydrophobic core, making it difficult for the TLP to leave the membrane. The three systems with the lowest barrier are combinations of TLP 135° PSC-AE with MP 180° PSC-NE (ΔΔGBT = 5.1 kJ/mol), TLP 135° PSC-NE with MP 180° PSC-NE (ΔΔGBT = 6.6 kJ/mol), and TLP 180° PSC-NE with MP 270° PSC-NE (ΔΔGBT = 7.2 kJ/mol).
Table 3

Free energy differences for important points during the translocation of TLP (adsorbed (A) and transmembrane (TM) states and profile maxima and minima) relative to free energy in solution and translocation free energy barriers, ΔΔGBT, for all simulated systems using the phenomenological model [in kJ/mol]

TLPMPΔGTMΔGAΔGMAXΔGMINΔΔGBT
90°none59.5NA59.60.059.6
PSC-NE360°63.6NA63.60.063.6
270°42.7NA42.70.042.7
180°22.5NA23.20.023.2
90°none58.9NA58.90.058.9
PSC-AE360°62.4NA62.40.062.4
270°40.4NA40.40.040.4
180°21.3NA21.30.021.3
135°none38.8NA39.40.039.4
PSC-NE360°41.9NA42.00.042.0
270°22.3NA24.10.024.1
180°4.5NA6.60.06.6
135°none37.4NA37.40.037.4
PSC-AE360°39.5NA39.50.039.5
270°20.4NA20.40.020.4
180°2.5−2.62.5−2.65.1
180°none17.1NA19.10.019.1
PSC-NE360°18.0NA20.00.020.0
270°1.5NA7.20.07.2
180°−13.5NA2.1−13.515.5
180°none14.6−6.514.6−6.521.1
PSC-AE360°15.7−5.415.7−5.421.1
270°−0.1−8.11.6−8.19.7
180°−15.5−16.61.3−16.617.8
270°none−31.6NA4.6−31.636.2
PSC-NE360°−31.7NA5.9−31.737.6
270°−46.8NA3.2−46.850.0
180°−46.7NA1.6−46.748.2
270°none−41.3−46.01.1−46.047.0
PSC-AE360°−39.4−44.91.1−44.946.0
270°−48.4NA1.0−48.749.8
180°−51.0NA0.8−51.352.1

The free energy of the adsorbed state is only shown when (meta)stable. The error was estimated to be below 2.5 kJ/mol, based on the profile asymmetry, and MPs are in increasing order of hydrophilicity.

Free energy differences for important points during the translocation of TLP (adsorbed (A) and transmembrane (TM) states and profile maxima and minima) relative to free energy in solution and translocation free energy barriers, ΔΔGBT, for all simulated systems using the phenomenological model [in kJ/mol] The free energy of the adsorbed state is only shown when (meta)stable. The error was estimated to be below 2.5 kJ/mol, based on the profile asymmetry, and MPs are in increasing order of hydrophilicity. Fig. 4 contains the interpolation of the translocation free energy barrier, ΔΔGBT, between the studied systems, separately for PSC-NE and PSC-AE TLPs. Both interpolation maps show a valley of low free barriers connecting a combination of TLP 135° × MP 180° in the lower part of the map with a combination of TLP 180° × MP 270° in the center of the map. This low barrier valley shows combinations of TLP and MP that are optimal for the translocation. Leaving this area by changing the hydrophobicity/hydrophilicity of either the TLP or the MP leads to an increase in the free energy barrier and more difficult translocation.
Figure 4

Interpolated translocation free energy barriers (ΔΔGBT) depending on translocating peptide (TLP) hydrophobic patch size (x axis) and transmembrane protein/peptide (MP) hydrophobic patch size (y axis). The translocation free energy barrier is the difference between the maximal and minimal free energy in the free energy profile. (A) and (B) show approximate ΔΔGBT for TLPs with fully hydrophilic termini (NE) and partially hydrophobic termini (AE), respectively. The free energy barrier for peptide translocation increases from purple to red. Black x-markers highlight MP × TLP combinations, which were simulated. Contour lines are drawn every 10 kJ/mol. To see this figure in color, go online.

Interpolated translocation free energy barriers (ΔΔGBT) depending on translocating peptide (TLP) hydrophobic patch size (x axis) and transmembrane protein/peptide (MP) hydrophobic patch size (y axis). The translocation free energy barrier is the difference between the maximal and minimal free energy in the free energy profile. (A) and (B) show approximate ΔΔGBT for TLPs with fully hydrophilic termini (NE) and partially hydrophobic termini (AE), respectively. The free energy barrier for peptide translocation increases from purple to red. Black x-markers highlight MP × TLP combinations, which were simulated. Contour lines are drawn every 10 kJ/mol. To see this figure in color, go online. Additionally, we have quantified the effect of MPs on membrane structure. We calculated the lipid defect for a pure membrane bilayer and for five membrane systems with various MPs (360, 270, 240, 210, 180° hydrophobicity) but no translocating peptide. The defect was defined as the average number of lipid headgroups in the membrane core located around the MP (see the Materials and methods for details). As expected, the lipid defect increased with increasing hydrophilicity of the MP present in the membrane (see Table S3), meaning that more hydrophilic MPs caused a larger membrane disruption. The MP with no hydrophilic patch (360°) actually caused a smaller lipid defect than in the pure membrane. As can be seen in Fig. S12, we found a strong negative correlation (r < −0.93) between lipid defects and ΔGTM values in all systems. The exceptions were systems with TLPs 270° (both PSC-NE and PSC-AE), for which the correlation was weaker (r = −0.70 and −0.83, respectively). Despite the small number of data points, these results suggest a connection between the ability of the MP to locally disrupt the membrane structure and the decrease in free energy for TLPs in the membrane. As anticipated, the membrane disruption seems to play a lesser role in the translocation of highly hydrophobic TLPs. Similarly to our simulations with the phenomenological model, we investigated the translocation of an amphiphilic α-helical peptide (TLP) using a Martini coarse-grained model (20, 21, 22). The TLP translocated across a phospholipid membrane along a symmetrical α-helical transmembrane protein/peptide (MP). All simulated proteins/peptides had a hydrophilic patch composed of serines and a hydrophobic patch composed of leucines. We studied seven MPs and three TLPs with various amphiphilicities, i.e., different sizes of their hydrophobic and hydrophilic patches. The individual sequences are listed in Tables 1 and 2. The peptide labels capture the hydrophilicity by the number of serine residues. TLPs LS5 and LS13 are polyleucines with 5 and 13 serine residues, respectively. The same labeling is applied to MPs, in which the extra three letters “SAG” at the beginning of the label correspond to the residues at the peptide/protein termini. SAGLS15 thus represents a polyleucine with SAG sequence at both termini and the total number of 15 serine residues. The calculated free energy profiles of TLP translocation are in Figs. S13–S15. We observed that with increasing hydrophilicity of the MP, the free energy of TLP in the membrane decreased. The studied MPs thus facilitated the insertion of the translocating peptide into the membrane and also stabilized its transmembrane state. This trend is clearly depicted in Fig. 5 with the translocation free energy profiles of a representative TLP in the presence of various MPs. TLP translocated along the MP following the same pathway as with our phenomenological model, i.e., from TLP adsorption on the membrane with an orientation parallel to the membrane plane, through the insertion of the peptide terminus into the membrane, and to the transmembrane state of TLP, where the TLP was oriented perpendicular to the membrane plane and with its hydrophilic patch oriented toward the MP (see Fig. S16).
Figure 5

Free energy profile of peptide LS9 translocating in the presence of various transmembrane proteins/peptides (MPs). MPs differed in their hydrophilic patch size: SAGLS5, SAGLS7, SAGLS11, and SAGLS15 had a total of 5, 7, 11, and 15 serine residues, respectively. Helical wheels of corresponding MPs are shown at the top. Green circles show hydrophilic residues, and hydrophobic residues are gray. ΔGIC and ΔGIN are the local free energy maxima of the C- and N-terminus insertion, respectively. ΔGTM is the free energy of TLP in the transmembrane state, and ΔGA is the free energy of the adsorbed state. The reference state for the above free energy differences is the peptide in solution. The magenta arrow shows the membrane free energy barrier (ΔΔGBM) of LS9 in the presence of SAGLS15. ΔΔGBM is defined as the difference between the local maximum and minimum in the membrane part of the free energy profile. The dark-cyan arrow shows the free energy barrier (ΔΔGBT) of the LS9 translocation across the membrane including the adsorption/desorption process. ΔΔGBT is defined as the difference between the global maximum and minimum in the free energy profile for peptide translocation. The error of profiles was estimated to be below 5 kJ/mol based on the difference between solution states (profile ends). To see this figure in color, go online.

Free energy profile of peptide LS9 translocating in the presence of various transmembrane proteins/peptides (MPs). MPs differed in their hydrophilic patch size: SAGLS5, SAGLS7, SAGLS11, and SAGLS15 had a total of 5, 7, 11, and 15 serine residues, respectively. Helical wheels of corresponding MPs are shown at the top. Green circles show hydrophilic residues, and hydrophobic residues are gray. ΔGIC and ΔGIN are the local free energy maxima of the C- and N-terminus insertion, respectively. ΔGTM is the free energy of TLP in the transmembrane state, and ΔGA is the free energy of the adsorbed state. The reference state for the above free energy differences is the peptide in solution. The magenta arrow shows the membrane free energy barrier (ΔΔGBM) of LS9 in the presence of SAGLS15. ΔΔGBM is defined as the difference between the local maximum and minimum in the membrane part of the free energy profile. The dark-cyan arrow shows the free energy barrier (ΔΔGBT) of the LS9 translocation across the membrane including the adsorption/desorption process. ΔΔGBT is defined as the difference between the global maximum and minimum in the free energy profile for peptide translocation. The error of profiles was estimated to be below 5 kJ/mol based on the difference between solution states (profile ends). To see this figure in color, go online. To describe the efficiency of peptide translocation, we defined a free energy barrier of translocation, ΔΔGBT, as the difference between the maximal and minimal free energy on the translocation profile, as in our phenomenological model. However, the translocation free energy barrier, ΔΔGBT, was not significantly affected by the presence of MP in most systems. As shown in Table 4, ΔΔGBT only decreased with increasing hydrophilicity of the MP for systems with very hydrophilic TLP (LS15) and very hydrophobic MPs (SAGL and SAGLS5). In contrast, the barrier, ΔΔGBT, increased with increasing hydrophilicity of the MP for systems with very hydrophobic TLP (LS5) and hydrophilic MPs (SAGLS7 to SAGLS15). In most systems, ΔΔGBT remained roughly constant. In these systems, ΔΔGBT corresponded to the difference between the TLP free energy in the solution and the adsorbed state, neither of which was affected by the presence of the MP.
Table 4

Free energy differences for important points on the TLP translocation pathway relative to the free energy in solution [kJ/ mol] for all systems simulated with the Martini model

TLPMPΔGAΔGICΔGTMΔGINΔΔGBTΔΔGBM
LS13none−104−2−1539143143
SAGL−106−13−2922129129
SAGLS5−104−27−4611116116
SAGLS7−106−43−68−5107101
SAGLS9−105−52−75−1410692
SAGLS11−105−57−83−2010986
SAGLS13−107−65−88−2610983
SAGLS15−106−62−85−2610881
LS9none−131−49−64−21132109
SAGL−134−61−75−32135103
SAGLS5−131−72−92−4413488
SAGLS7−131−86−113−5613375
SAGLS9−133−93−120−6913364
SAGLS11−131−96−126−7113462
SAGLS13−133−103−133−8013454
SAGLS15−135−106−132−8313552
LS5none−148−95−124−8714861
SAGL−149−99−127−9114958
SAGLS5−149−108−142−10015050
SAGLS7−150−117−154−10815446
SAGLS9−147−121−157−11316043
SAGLS11−149−126−162−11516347
SAGLS13−151−130−165−12216644
SAGLS15−150−130−164−12016444

The error was estimated to be below 5 kJ/mol based on the profile asymmetry. MPs are listed in increasing order of hydrophilicity.

Free energy differences for important points on the TLP translocation pathway relative to the free energy in solution [kJ/ mol] for all systems simulated with the Martini model The error was estimated to be below 5 kJ/mol based on the profile asymmetry. MPs are listed in increasing order of hydrophilicity. To capture the effect of the MP on the translocation of TLPs, which were strongly adsorbed on the membrane, we defined the membrane free energy barrier, ΔΔGBM. ΔΔGBM was calculated as the difference between the maximal and minimal free energy of the TLP in the membrane, i.e., the membrane insertion part of the translocation profile. ΔΔGBM decreased linearly with increasing hydrophilicity of the MP until a threshold was reached. Any subsequent increase in the MP hydrophilicity did not lead to a further decrease of the barrier height. For visualization of the results, we also constructed an interpolated surface/map of the ΔΔGBM (see Fig. 6). The largest decrease was for the most hydrophilic MP (SAGLS15), which decreased the membrane barrier to roughly half for the MPs LS9 and LS13 compared with no MP. The decrease was “only” about one-third for the most hydrophobic TLP we investigated (LS5).
Figure 6

Interpolated height of membrane free energy barrier (ΔΔGBM) depending on translocating peptide (TLP) hydrophobicity (x axis) and transmembrane protein (MP) hydrophobicity (y axis). The employed hydrophobicity scale is based on the partitioning of N-acetyl-amino-acid amides in octanol (see the Materials and methods for details). Black x-markers indicate simulated MP × TLP combinations. Contour lines are drawn every 10 kJ/mol. To see this figure in color, go online.

Interpolated height of membrane free energy barrier (ΔΔGBM) depending on translocating peptide (TLP) hydrophobicity (x axis) and transmembrane protein (MP) hydrophobicity (y axis). The employed hydrophobicity scale is based on the partitioning of N-acetyl-amino-acid amides in octanol (see the Materials and methods for details). Black x-markers indicate simulated MP × TLP combinations. Contour lines are drawn every 10 kJ/mol. To see this figure in color, go online. To characterize the effect of MPs on a membrane, we calculated the water defect in the lipid bilayer around each simulated MP without the presence of TLP. The water defect was defined as the average number of water beads in a cylinder around MP (see the Materials and methods for details). All simulated MPs caused a larger and deeper water defect than was observed in the pure membrane (see Table S4). The water defect increased with increasing hydrophilic patch size of MP. The obtained values of water defect almost perfectly negatively correlated (r = −0.99) with the hydrophobicity of the MP. We also identified a high negative correlation of water defects with values of ΔGIC, ΔGIN, and ΔGTM (r < −0.94) for all three simulated TLPs (see Fig. S17, A–C). We also observed a high negative correlation between water defects and ΔΔGBM (r < −0.90) (see Fig. S17 D). All these correlations suggest that the MP’s disruption of the membrane structure could have a direct impact on the translocation free energy, with more hydrophilic MPs causing a larger effect. Note that no water pores (continuous water channels) were detected around any of the MPs during our simulations. To test whether there is a causal relation between the membrane disruption and the ability of translocating peptide to insert itself into the membrane, we simulated membranes with three different artificially created water defects (see the Materials and methods for details). The artificial water defects closely corresponded to the water defects observed around MPs SAGL, SAGLS7, and SAGLS15. The calculated free energy profiles of LS9 peptide translocating through a membrane with a defect and through a membrane with the corresponding MP were very similar (see Fig. S18).

Discussion

We investigated the effect of transmembrane proteins/peptides (MPs) with various properties on the translocation of amphiphilic peptides (TLPs) across the membrane using 1) a phenomenological model of amphiphilic peptides and 2) a Martini model with serine-leucine peptides. We focused on the translocation of a single peptide across the membrane because we were interested in translocation mechanisms at a low peptide concentration. In both employed models, MPs with a hydrophilic patch caused a decrease in the free energy of TLP in the membrane. In other words, TLPs inserted themselves into the bilayer more easily and were also more stable in the transmembrane state in the presence of such MPs. Increasing the hydrophilic patch size of the MP led to stronger stabilization of the TLP in the membrane, both during the insertion and in the transmembrane state. In some cases, the stabilization of the TLP was so strong that the TLP preferred its transmembrane state to its adsorbed state. The pathway of peptide translocation in the presence of an MP was the same as the translocation through the lipid membrane without an MP, as described in previous studies (18,36, 37, 38). The TLP first adsorbed itself onto the membrane in an orientation parallel to the membrane plane. As the TLP inserted into the membrane, it changed its orientation. In the transmembrane state, the peptide was positioned in an orientation perpendicular to the membrane plane. During the translocation, the hydrophilic patches of the TLP and MP were usually facing each other. This interaction presumably stabilized the peptide during the insertion stage and in the transmembrane state. Usually, by increasing the hydrophilicity of the MP, TLP translocation became easier as the free energy barrier of either the entire translocation profile or at least the membrane part of the profile decreased. In our Martini simulations, the limiting step of the TLP translocation was usually the process of desorption from the membrane, and the TLP was usually more stable in the membrane than in the solvent. However, note that the free energy of peptide desorption from the membrane was likely overestimated because of the lack of unfolding free energy of the peptide in solution, which cannot be calculated using either of the employed models. We can only estimate the unfolding free energy using a generic sequence-unspecific contribution of 1.7 kJ/mol per residue, leading to a free energy decrease of roughly 35 kJ/mol (favoring solution state) (39). Another significant barrier in the peptide translocation profile was connected to the partial insertion of the peptide into the membrane. Although the adsorption/desorption process was not affected by the presence of the studied MPs, membrane insertion barriers were modulated by it. In our Martini simulations, the free energy barrier for crossing the membrane always decreased in systems containing any of the simulated MPs compared with systems with the TLP passing through the membrane alone. This means that every simulated MP enhanced the process of the TLP crossing the membrane. In most of our simulations with phenomenological model, the limiting step of the translocation was the TLP passing through the lipid bilayer. For specific MP × TLP combinations, the TLP preferred the membrane core environment to the solution. In such cases, the loss of the membrane insertion barrier caused the process of TLP leaving the membrane core to be the only limiting step of the translocation. We refer to this situation as overstabilization of the TLP in the membrane. The comparison of results from the phenomenological model and Martini force field is not straightforward, because the phenomenological model does not correspond to any specific sequence. Nevertheless, if we wanted to align the interpolated barrier maps (Fig. 4 A; Fig. 6), TLP LS5 would have to correspond to PSC-NE with a hydrophobic patch size of 180°, and MP SAGLS5 would have to correspond to PSC-NE with a much larger hydrophobic patch. However, LS5 is more hydrophobic than SAGLS5, making this requirement contradictory. Therefore, the interpolated barrier maps are not comparable between the models. Based on our analysis of membrane defects (early stage of pore formation), we argue that the discrepancy between the employed models mostly stems from the different pore-forming abilities of the simulated lipid bilayers. It is known that Cooke-Deserno lipid bilayers, which we used in our phenomenological simulations, overestimate the rate of lipid flip-flop by several orders of magnitude (40). This causes an increased formation of short-lived pores in the bilayer, especially near the embedded MP, which was able to stabilize lipid heads in the pore. In contrast, the Martini model is known to overestimate the free energy of pore formation compared with atomic simulations (41). Therefore, it is expected that the same hydrophilicity of the MP caused larger structural changes in the membrane described by the Cooke-Deserno model and hence influenced the translocation more than in the Martini model. However, note that the direct relation between the membrane defect and easier peptide translocation could be more complex for specific sequences. For example, the presence of negatively charged residues in the MP is expected to affect the translocation of positively and negatively charged peptides differently as a result of electrostatic interactions. The occurrence of membrane defects around the simulated MPs suggests that MPs with hydrophilic patches have properties similar to scramblases: proteins enabling the spontaneous translocation (scrambling) of phospholipids between the leaflets of the membrane. The existence of a hydrophilic pathway on the surface of a transmembrane segment seems to be critical for lipid scrambling, as seen with TMEM16, which contains a groove lined with hydrophilic residues that lipids translocate along (10), and opsin, whose hydrophilic pathway is formed by a cooperation between two α-helices (42). Based on our findings, we propose that MPs that are able to scramble lipids can also act as translocation enhancers facilitating the translocation of amphiphilic peptides. AMPs or CPPs could then use naturally occurring scramblases or proteins with similar structural features to enter cells more effectively. In the future, the specific targeting of potential bacterial scramblases could also be employed when designing new AMPs and CPPs.

Conclusions

We investigated the effect of transmembrane proteins/peptides with a hydrophilic patch on the translocation of amphiphilic peptides across a phospholipid membrane. Using coarse-grained Martini and phenomenological models, we consistently showed that the transmembrane proteins/peptides could enhance the translocation of peptides, including AMPs and CPPs. Increasing the hydrophilicity of the transmembrane proteins/peptides leads to a larger decrease in the insertion barrier and stronger stabilization of the translocating peptide in the membrane. Once the transmembrane protein/peptide was too hydrophilic, the translocating peptide became more stable in the transmembrane state than in the adsorbed state. In the phenomenological model, we even investigated transmembrane proteins/peptides, which hindered the translocation by overstabilizing the translocating peptide in the membrane. Nevertheless, our results demonstrate that the translocation of amphiphilic peptides across membranes could be enhanced by the presence of native membrane proteins that are not fully hydrophobic in their transmembrane part, such as scramblases. These findings also suggest that a mixture of two peptides could be more effective in transmembrane transport, when a very hydrophobic peptide with a small hydrophilic patch inserts into the membrane and then enhances the translocation of another more hydrophilic peptide.

Author contributions

L.B. carried out all simulations and analyzed the data. I.K. devised technical solutions to problems in simulations. R.V. designed the research. L.B., I.K., and R.V. wrote the article.
  31 in total

1.  Efficient, multiple-range random walk algorithm to calculate the density of states.

Authors:  F Wang; D P Landau
Journal:  Phys Rev Lett       Date:  2001-03-05       Impact factor: 9.161

2.  Canonical sampling through velocity rescaling.

Authors:  Giovanni Bussi; Davide Donadio; Michele Parrinello
Journal:  J Chem Phys       Date:  2007-01-07       Impact factor: 3.488

Review 3.  Cell-Penetrating Peptides: From Basic Research to Clinics.

Authors:  Giulia Guidotti; Liliana Brambilla; Daniela Rossi
Journal:  Trends Pharmacol Sci       Date:  2017-02-14       Impact factor: 14.819

Review 4.  Antimicrobial peptides: Promising alternatives in the post feeding antibiotic era.

Authors:  Jiajun Wang; Xiujing Dou; Jing Song; Yinfeng Lyu; Xin Zhu; Lin Xu; Weizhong Li; Anshan Shan
Journal:  Med Res Rev       Date:  2018-10-24       Impact factor: 12.944

5.  Optimal Hydrophobicity and Reorientation of Amphiphilic Peptides Translocating through Membrane.

Authors:  Ivo Kabelka; Robert Vácha
Journal:  Biophys J       Date:  2018-08-18       Impact factor: 4.033

6.  Binding and reorientation of melittin in a POPC bilayer: computer simulations.

Authors:  Sheeba J Irudayam; Max L Berkowitz
Journal:  Biochim Biophys Acta       Date:  2012-08-02

Review 7.  Cell-penetrating peptides: 20 years later, where do we stand?

Authors:  Chérine Bechara; Sandrine Sagan
Journal:  FEBS Lett       Date:  2013-05-10       Impact factor: 4.124

8.  Out-of-the-groove transport of lipids by TMEM16 and GPCR scramblases.

Authors:  Mattia Malvezzi; Kiran K Andra; Kalpana Pandey; Byoung-Cheol Lee; Maria E Falzone; Ashley Brown; Rabia Iqbal; Anant K Menon; Alessio Accardi
Journal:  Proc Natl Acad Sci U S A       Date:  2018-06-20       Impact factor: 11.205

9.  The energetics of transmembrane helix insertion into a lipid bilayer.

Authors:  Alan Chetwynd; Chze Ling Wee; Benjamin A Hall; Mark S P Sansom
Journal:  Biophys J       Date:  2010-10-20       Impact factor: 4.033

Review 10.  Membrane Active Peptides and Their Biophysical Characterization.

Authors:  Fatma Gizem Avci; Berna Sariyar Akbulut; Elif Ozkirimli
Journal:  Biomolecules       Date:  2018-08-22
View more

北京卡尤迪生物科技股份有限公司 © 2022-2023.