# Understanding Passive Membrane Permeation of Peptides: Physical Models and Sampling Methods Compared

Liuba Mazzanti, Tâp Ha-Duong

► To cite this version:

Liuba Mazzanti, Tâp Ha-Duong. Understanding Passive Membrane Permeation of Peptides: Physical Models and Sampling Methods Compared. International Journal of Molecular Sciences, 2023, 24 (5), pp.5021. 10.3390/ijms24055021. hal-05375544

HAL Id: hal-05375544

https://hal.science/hal-05375544v1

Submitted on 20 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

# Understanding Passive Membrane Permeation of Peptides: Physical Models and Sampling Methods Compared

Liuba Mazzanti and Tâp Ha-Duong*

BioCIS, CNRS, Université Paris-Saclay, 17 Avenue des Sciences, 91400 Orsay, France

* Correspondence: tap.ha-duong@universite-paris-saclay.fr

**Abstract:** The early characterization of drug membrane permeability is an important step in pharmaceutical developments to limit possible late failures in preclinical studies. This is particularly crucial for therapeutic peptides whose size generally prevents them from passively entering cells. However, a sequence-structure-dynamics-permeability relationship for peptides still needs further insight to help efficient therapeutic peptide design. In this perspective, we conducted here a computational study for estimating the permeability coefficient of a benchmark peptide by considering and comparing two different physical models: on the one hand, the inhomogeneous solubility-diffusion model, which requires umbrella-sampling simulations, and on the other hand, a chemical kinetics model which necessitates multiple unconstrained simulations. Notably, we assessed the accuracy of the two approaches in relation to their computational cost.

**Keywords:** peptide membrane permeability; molecular dynamics simulation; umbrella sampling; Markov State Model; free energy profile

## 1. Introduction

To reach an intracellular target, a drug must cross the cellular membrane. This process can be performed by endocytosis, which is generally involved in internalizing macro-molecules or drug nanocarriers [1], or by using membrane transporter proteins, such as the human ATP binding cassette (ABC) [2] and solute carrier (SLC) [3] families, or the outer membrane porins [4] and TonB-dependent transporters [5] of gram-negative bacteria. However, the predominant mechanism for xenobiotics is a passive diffusion across the lipid bilayer along a concentration gradient [6,7]. It is, thus, essential in pharmaceutical developments to characterize the drug membrane permeability to anticipate their pharmacokinetic properties. This is particularly crucial for peptide-based therapeutics. Indeed, peptides are very promising compounds for modulating protein-protein interactions (PPIs), especially those involved in intracellular signaling pathways [8]. This is partly due to their size, which is well-appropriate for covering the generally large molecular surfaces involved in PPIs. However, this size may turn out to be a disadvantage for membrane permeation.

Passive membrane permeability of drugs can be characterized experimentally by using the parallel artificial membrane permeability assay (PAMPA), which consists in measuring the number of compounds that have crossed a planar artificial membrane from donor to acceptor wells [9]. Alternative methods have been developed to monitor in real-time the entry of drugs into unilamellar liposomes using fluorescent probes [10,11]. In both approaches, the drug membrane permeability *P* is quantified as *P* = *J*/*ΔC*, where *J* is the flux of the drug through the membrane and *ΔC* the difference between the donor (or outside) and acceptor (or inside) drug concentrations.

Throughout a membrane permeation process, a peptide can bind the lipid headgroup-water interface, fold into secondary structures, form aggregates, insert into the lipid hydrophobic tails, and perturb the bilayer organization [12,13]. These multiple possible events depend on the peptide apolar/polar/charged amino acid composition, its structure and

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Mazzanti, L.; Ha-Duong, T. Understanding Passive Membrane Permeation of Peptides: Physical Models and Sampling Methods Compared. *Int. J. Mol. Sci.* **2023**, *24*, 5021. https://doi.org/10.3390/ijms24055021

Academic Editors: Geoffrey Brown, Enikö Kallay and Andrzej Kutner

Received: 1 February 2023

Revised: 21 February 2023

Accepted: 23 February 2023

Published: 6 March 2023

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
flexibility, the solvent pH, the ionic strength, and the membrane content in lipids. All these factors make it challenging to predict the membrane permeability and deeply understand the mechanism of peptide translocation. To tackle this issue, molecular dynamics (MD) simulations of peptide crossing lipid bilayers are valuable tools to gain insight into the physical-chemical factors that govern the permeation process.

However, membrane permeation occurs on a millisecond to second timescale [10] which are generally out of reach for classical MD simulations. Thus, enhanced sampling methods, such as metadynamics or umbrella sampling (US), are needed to provide valuable information about the compound propensity to partition between water and membrane, the structures of the intermediate states within the bilayer, and about the mechanisms of permeation [14,15]. In particular, trajectories from US simulations allow to compute drug free energy profiles (FEP) across lipid bilayers and can be employed to estimate their membrane permeability by using the so-called inhomogeneous solubility-diffusion model (ISDM) proposed by Marrink and Berendsen [16]. This approach has been applied to several drugs over the past few years [14,17–19] and recently to peptides [20].

Since membrane permeability is a kinetic property, it can be naturally computed by using a kinetic model of the permeation process. This consists in describing the elementary steps of a compound membrane permeation and formulating each of the corresponding reaction rates as a function of the concentrations in the intermediate states. Then, the master differential equations are solved to yield the compound concentrations on both sides of the membrane as a function of time, from which the membrane permeability can be retrieved with Fick's law of diffusion [18]. Importantly, this method requires to preliminary determine the rate constants of the kinetic model. This can be done with Markov State Model (MSM) analyses which can extract long-time kinetic information of a molecular system from an ensemble of "short" MD trajectories starting from multiple initial configurations [21,22].

For small compounds, calculations based on kinetic models appear to provide mem-
brane permeability values that better correlate with experiments than those based on
ISDM [18]. Nevertheless, as far as we know, the former approach has never been applied to
peptides. In the present study, we comparatively assessed the two computational methods
to estimate the membrane permeability for a benchmark cyclic peptide of 10 amino acid
residues. Notably, we scrutinized the accuracy of the two approaches in relation to their
computational cost.

## 2. Results and Discussion

We applied the two computational approaches to a paradigmatic ten amino acids cyclic peptide for which PAMPA assays were performed [9]. Constrained or unconstrained MD trajectories of the peptide assembled with a planar lipid bilayer were collected to ultimately calculate the peptide permeability coefficient through Equation (4) or (8) for the ISDM or the MSM-based model, respectively. In both approaches, the free energy profile of the peptide perpendicularly crossing the planar bilayer membrane is crucial for estimating its permeability coefficient. In the ISDM method, the latter is directly related to the FEP ΔG(z) through Equations (2)–(4). On the other hand, in the MSM-based approach, without the knowledge of the energy barrier location, the Markov State Models could be biased since the phase space sampling by short trajectories depends on the starting points along the collective variable. Therefore, the initial configurations must correspond to maximum energy states for the trajectories not to be trapped in energy minima. Thus, we will present first the free energy analyses and, subsequently, the quantitative estimates for the permeability coefficients.

2.1. *Free Energy Profiles*

The weighted histogram analysis method (WHAM) [23] applied to the umbrella
sampling trajectories for the cyclic decapeptide designated as CDP5 (Section 3.1) yields the
free energy profile displayed in Figure 1. Moving from the bulk water to the membrane, the
peptide encounters a first energy barrier of about $\Delta G_{in} = 3$ kcal/mol located on the surface of the lipid headgroups in the water phase ($z = 2.4$ nm). Then the FEP exhibits a minimum on the other side of the lipid headgroup plane, inside the membrane ($z = 1.4$ nm), where the amphiphilic peptide is stabilized by both hydrophobic and polar interactions with the lipid tails and headgroups, respectively. We note that the free energy of this stable state is roughly equal to that one in the bulk water. Therefore, the free energy barrier for moving out from the interior of the membrane into the bulk water is also $\Delta G_{out} = \Delta G_{in} = 3$ kcal/mol. The major energy barrier for crossing the membrane is located at its center ($z = 0.0$ nm), indicating that the rate-limiting step during the peptide permeation is the flip-flop passage from one lipid layer to the other. The necessary energy to overcome this barrier is about $\Delta G_{flip} = 12$ kcal/mol.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Symmetrized CDP5 free energy profile computed from umbrella sampling trajectories. Error bars were estimated using the gmx wham bootstrap algorithm [24]. The orange vertical line and shaded area indicate the mean and one $\sigma$ of the lipid headgroup positions, respectively. Free energy difference $\Delta G(z) = G(z) - G(z_{max})$ was calculated with $z_{max} = 3.8$ nm in the water phase. Four representative structures of the peptide-membrane system are shown for the peptide positions $z = 0.0, 1.4, 2.0$, and 2.7 nm. Peptides, lipid headgroups, and lipid tails are displayed using sticks, semi-transparent spheres, and lines, respectively.</figcaption></figure>

The FEP computed for the CDP5 peptide is very similar to those of amphiphilic cyclic hexapeptides reported by Sugita et al. [20], indicating that these cyclic peptides of comparable size should have the same mechanism of membrane translocation in which the flip-flop event within the bilayer is the rate-limiting step. This contrasts with the FEPs of amphiphilic small compounds computed by Dickson et al. [18], which exhibit two marked energy minima located at the inner lipid headgroup surfaces much lower than the free energy in the bulk water. For these small compounds, the energy barrier at the membrane center is generally small, depending on their hydrophobicity, indicating that the major rate-limiting step of their permeation process is moving from the interior of the membrane into the bulk water.

To compute the peptide membrane permeability using kinetic models, we built a first Markov State Model (CDP5_2.4) from unconstrained MD simulations of the peptide with initial positions at $z = 2.4$ nm and $z = 0.0$ nm corresponding to the two maxima of its FEP computed from US trajectories. Moreover, to assess the sensitivity of this approach to the initial configurations, we also built Markov State Models (CDP5_2.3 and CDP5_2.6), respectively from trajectories starting at $z = 2.3$ nm and $z = 0.0$ nm on the one hand, and
from trajectories starting from *z* = 2.6 nm and *z* = 0.0 nm, on the other hand. A fourth Markov State Model (CDP5_all) was also built using all the aforementioned trajectories (Section 3.3). From each of these sets of unconstrained MD trajectories, we computed the peptide stationary distributions *π* and the free energy profiles *ΔG* = −*k*<sub>B</sub>*T* log *π*, which are shown in Figure 2.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** CDP5 free energy profiles computed from four different sets of unconstrained trajectories (Section 3.3). Diamonds indicating the discrete states are colored according to the metastable state to which they belong. FEPs were set to zero in the water phase, including the FEP computed from US trajectories (black lines).</figcaption></figure>

FEPs derived from unconstrained trajectories show significant differences with
respect to the profile computed from the US trajectories: First, the energy barriers
to cross the lipid headgroups are smaller (by about 2 kcal/mol) than the US ones,
except for the (CDP5_2.6) set of trajectories. Secondly, the positions of the energy
minima inside the membrane are shifted toward the lipid tails (z = 1.0 nm) compared
to the US FEP (z = 1.4 nm). However, the most striking difference is the free energy
barrier corresponding to the flip-flop step: the unconstrained simulations yield a central
free energy barrier of about 2 kcal/mol, while the US flip-flop barrier almost reaches
12 kcal/mol. This suggests that the peptide membrane permeability estimated with the
kinetic models will be much higher than in the ISDM-based approach.

It is interesting to note that the four FEPs generated by the four sets of unconstrained
simulations display different qualitative behaviors, implying that a slight change in their
initial conditions could entail a sensible variation in the results. In particular, we notice that
the free energy needed to translocate the peptide from inside the membrane into the water
phase is larger than the energy barrier associated with the flip-flop step in the CDP5_2.6
model, while it is smaller in the three other cases (Figure 2). The emergence of high in-out
free energy barriers in the CDP5_2.6 case is a consequence of the initial peptide positions,
which were slightly shifted from the minor maxima of the US FEP (z = 2.4 nm) toward the
water bulk. Therefore, fewer peptide trajectories could naturally sample this maximum
energy state which increases the barrier height, compared to the three other cases (where it
could even be observed a local shallow energy minimum at this location). This highlights
the sensitivity of the peptide FEPs computed from unconstrained MD simulations upon
their initial configurations.
## 2.2. Membrane Permeability Calculations

In the framework of the ISDM approach, an estimate of the peptide permeability coefficient can be directly obtained from the US trajectories by extracting the integrated autocorrelation times $\tau_{zz}$, the collective variable variance $\sigma_z^2$, and the free energy profile $\Delta G(z)$ and applying Equations (2)–(4). For the CDP5 peptide, we get the value reported in Table 1. Compared to experiments, the ISDM-based estimation of the peptide membrane permeability is about one order of magnitude lower than the value found in PAMPA assays [25].

<table><caption>**Table 1.** Membrane permeability coefficient results for the cyclic decapeptide CDP5. PAMPA value was taken from [25].</caption><thead><tr><th>Method</th><th>log[P(cm/s)] ± σ<sub>P</sub></th></tr></thead><tbody><tr><td>PAMPA</td><td>-6.7</td></tr><tr><td>ISDM</td><td>-8.1 ± 0.4</td></tr><tr><td>MSM CDP5_2.3</td><td>-1.0 ± 0.1</td></tr><tr><td>MSM CDP5_2.4</td><td>-1.2 ± 0.1</td></tr><tr><td>MSM CDP5_2.6</td><td>-1.4 ± 0.1</td></tr><tr><td>MSM CDP5_all</td><td>-1.0 ± 0.1</td></tr></tbody></table>

On the other hand, further calculations and analysis are necessary to get an estimate of the permeability coefficient within the liposome kinetic model approach. Rate constants are derived as the inverse of the mean first passage times (MFPTs) and then converted to apparent rate constants for the liposome model using Equation (6). As detailed in Appendix A, integration of the kinetic model Equation (5) yields a three-exponential time evolution for the substrate in the inner aqueous compartment (Equation (A7)). By substituting $k_I$ of Equation (8) with the fastest rate constants $\kappa_n$ appearing in the substrate time evolution, we obtained the permeability coefficients reported in Table 1. Compared to the PAMPA value, it appears that the liposome kinetic models overestimate the peptide membrane permeability by about five to six orders of magnitude.

While comparing theoretical with experimental estimates of the peptide membrane permeability, we should bear in mind that PAMPA assays use hexadecane plates separating two water compartments [25] instead of a lipid bilayer or spherical liposomes. Thus, if we assume that crossing a hexadecane plate is easier than a lipid bilayer, this difference in the membrane nature might explain the lower permeability found by the ISDM-based method compared to experiments. On the other hand, the permeability estimated with liposome kinetic models being proportional to the liposome radius, it should be intuitively much lower than the permeability through a planar membrane. This is the opposite tendency that we can notice in Table 1 and, altogether, it appears that the MSM-based approach greatly overestimates the peptide membrane permeability compared to experiments.

The discrepancy between the ISDM-based and the MSM-based approaches can be understood as a consequence of the difference in the free energy barriers of the FEPs computed by the two methods. As the energy barrier for the flip-flop step is larger in the US FEP than in the unconstrained simulations-derived one, permeation through the membrane is intuitively less easily achieved within the ISDM framework than from the standpoint of the MSM-based liposome kinetic models. To a lesser extent, the slight variations of the energy barrier between the MSM-based FEPs can account for the slight differences in the permeability coefficient estimations (Table 1).

# 3. Methods and Materials

## 3.1. *Building the Peptide-Membrane System*

The peptide that we chose for the benchmark of the two methods that compute membrane permeability is a cyclic decapeptide named CDP5 (Figure 3). Four out of its ten residues were N-methylated to improve its cell permeability and oral bioavailability [26,27]. Its membrane
permeability coefficient was measured using PAMPA experiments [9] with a hexadecane layer separating the donor and acceptor compartments [25]. It should be noted that this hexadecane layer differs from the lipid bilayer that will be considered in our simulations and which is a more realistic model of experimental liposomes. Nevertheless, as far as this study is concerned, the experimental value of CDP5 permeability coefficient reported in [25] (log P<sub>exp</sub> = -6.7) will be used as a reference for our theoretical results.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Sequence and structure of the cyclic decapeptide CDP5 with four N-methylated residues.</figcaption></figure>

The modified peptide CDP5 was parameterized using the highly optimized chemical building blocks from CGenFF [28] to be compatible with the CHARMM36 force field for lipids [29]. An initial three-dimensional structure was generated from its sequence using OpenBabel [30] and PyMol [31]. After an energy minimization of 5000 steps and two short equilibration runs (1 ns in NVT and 2.5 ns in NPT ensemble), CDP5 was submitted to a 100 ns MD simulation in water with a sodium chloride concentration of 0.15 mol/L, using a V-rescale thermostat and a Parrinello-Raman barostat set at 300 K and 1 bar, respectively. Lennard-Jones potentials were cut off at 1.2 nm, electrostatic interactions were calculated using the smooth PME method [32], and covalent bonds with hydrogen were constrained using LINCS algorithm [33], for all simulations in this work. The resulting trajectory showed strong structural stability of the CDP5 peptide and was clustered based on RMSD with a cutoff value of 0.15 nm. The representative structure of the most populated cluster (Figure 3) provided the peptide conformation to be assembled with the membrane.

A lipid bilayer composed of 25 POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) molecules per leaflet were built at the temperature of T = 293.15 K using CHARMM-GUI [34]. The membrane was energy minimized and equilibrated in water with 0.15 mol/L of NaCl by running a 250 ps NVT and a 125 ps NPT equilibration with a semi-isotropic pressure scaling, followed by three 500 ps NPT runs with decreasing position restraints on lipid phosphate groups from 1000 to 0 kJ/mol/Å². Then, the peptide and the lipid bilayer were assembled in a rectangular box, solvated, energy minimized, and equilibrated with 1 ns NVT and 1 ns NPT runs by decreasing the position restraints on peptide backbone heavy atoms and lipid phosphate groups from 4000 and 1000 to 50 and 0 kJ/mol/Å², respectively. The final conformation of the equilibration procedure was taken as the initial configuration for US simulations.

GROMACS 2021 [35] was used for all simulations, with CHARMM36 [29] force-field for lipids, and TIP3P [36] model for water.

## 3.2. ISDM-Based Calculation of Membrane Permeability

This subsection presents the calculations that were performed to estimate the peptide membrane permeability by using the inhomogeneous solubility-diffusion model (ISDM) [16].

### Umbrella Sampling and WHAM Analysis

In the initial configuration of the peptide-membrane system, the solute is located in the bulk water at a distance z<sub>max</sub> = 4 nm from the center of the lipid bilayer (the membrane plane
being oriented perpendicular to the z-axis). Starting from this initial position, the solute was pulled toward and inside the lipid membrane by an increment of Δz = 0.1 nm until it reached the position z_min = 0 nm corresponding to the center of the bilayer. For each peptide position (or window), which is maintained with an umbrella potential of 1000 kJ/mol/Å², energy minimization was performed, followed by a 1 ns NVT and 1 ns NPT equilibration. Then, a production simulation of 250 ns was performed. The final conformation of the equilibration stage was pulled and used as a starting conformation for the next window.

The last 130 ns of the production trajectories were analyzed with the WHAM analysis [23] implemented in GROMACS [24]. More specifically, taking as input the values of the collective variable z for each window, gmx wham tool was used to compute both the free energy profile ΔG(z) = G(z) - G(z_max) and the normalized integrated autocorrelation time τ_zz defined as:

$$ \tau_{zz} = \frac{1}{\sigma_z^2} \int_0^\infty [z(0) - \langle z \rangle][z(t) - \langle z \rangle]dt \quad (1) $$

where σ_z² is the variance of the collective variable z defined as σ_z² = ⟨z²⟩ - ⟨z⟩². It could be noted that gmx wham provides both averages and standard errors of the FEP and τ_zz by using bootstrap calculations [24].

## Inhomogeneous Solubility-Diffusion Model

The ISDM considers that, due to the inhomogeneous nature of lipid membranes, the diffusion rate of a solute strongly depends on its position in the latter [16]. Moreover, it assumes that the diffusion process obeys the classical diffusion equation, implying that the dependence of ⟨z²⟩ on time is linear. Thus, the position-dependent diffusion coefficient can be calculated as the ratio of the variance σ_z² over the integrated autocorrelation time τ_zz provided by WHAM analyses of US simulations:

$$ D(z) = \frac{\sigma_z^2}{\tau_{zz}} \quad (2) $$

Using the diffusion coefficient D(z) and the free energy ΔG(z), the local resistance of the solute permeation through the membrane can be defined as [16]:

$$ R(z) = \frac{\exp[\beta\Delta G(z)]}{D(z)} \quad (3) $$

where β = 1/(k_B T), with T the temperature and k_B the Boltzmann's constant. The sum of all the R(z) along the variable z gives the global resistance and its inverse yields the membrane permeability:

$$ P_{ISDM} = \left[ \int_{-z_{max}}^{z_{max}} R(z) dz \right]^{-1} \quad (4) $$

To estimate the errors on the permeability coefficient, 100 bootstrapped trajectories were generated from the original trajectory of each window using an in-house script. Then, for each bootstrapped trajectory, the permeability coefficient was calculated as described above. The final results were computed as the mean of these 100 values and errors were given by the standard deviation.

## 3.3. MSM-Based Calculation of Membrane Permeability

In this subsection, we present the different steps of the method for computing the peptide membrane permeability by using a kinetic Markov State Model (MSM) [21].

### Multiple Unconstrained MD Simulations

A kinetic Markov State Model is essentially based on the transition probabilities between the model discretized states. These probabilities can be estimated by running multiple
unbiased MD simulations, which make the system pass through its different minimum energy states. Nevertheless, to sufficiently sample the transition events, it is recommended to run these simulations from out-of-equilibrium initial configurations, which can "go down the hill" to the nearest metastable states in a short time [21,22]. The identification of these maximum energy states implies the preliminary determination of the system's free energy landscape. Regarding the permeation of CDP5 peptide, its FEP exhibits two main transition states, a major one at the center of the bilayer ($z = 0$ nm) and a minor one at the surface of the lipid headgroups (around $z = 2.4$ nm) (Figure 1). Consequently, we decided to run simulations of the peptide-membrane system, starting from four sets of 75 conformations each, extracted from the previous US simulations, with peptide positions $z = 0.0, 2.3, 2.4$, and $2.6$ nm. Each of these $4 \times 75$ initial configurations were submitted to a 100 ns unconstrained MD production in the same conditions as the US simulations.

## Building the Markov State Models

To study the sensitivity of the computed permeability to the initial configurations of the peptide-membrane system, we built four Markov models from four different sets of MD trajectories (Table 2).

<table><caption>**Table 2.** Markov State Models built from four different sets of 100 ns unconstrained trajectories of the peptide-membrane system.</caption><thead><tr><td>MSM Name</td><td>Nb of Trajectories</td><td>Initial Positions</td></tr></thead><tbody><tr><td>CDP5_2.3</td><td>2 × 75</td><td>$z = 0.0$ and $2.3$ nm</td></tr><tr><td>CDP5_2.4</td><td>2 × 75</td><td>$z = 0.0$ and $2.4$ nm</td></tr><tr><td>CDP5_2.6</td><td>2 × 75</td><td>$z = 0.0$ and $2.6$ nm</td></tr><tr><td>CDP5_all</td><td>4 × 75</td><td>$z = 0.0, 2.3, 2.4$, and $2.6$ nm</td></tr></tbody></table>

PyEMMA [22] was used to build and analyze the four MSMs. For each kinetic model, the estimated implied relaxation timescale was plotted as a function of different lag times $τ$ (Figure 4). From these graphs, we considered that the implied relaxation timescales are approximately constant from 3 timesteps, i.e., $τ = 300$ ps, and that this lag time is suitable for calculating accurate kinetic coefficients of all MSMs. These plots also indicate that the models are characterized by three implied slow relaxation timescales (blue, red, and green lines in Figure 4). The transition matrix was estimated using the above lag time value, and a Chapman-Kolmogorov test was performed to verify the Markovianity of each model [22].

Next, PCCA+ [37] clustering was applied to describe all trajectories of each MSM in terms of a few metastable states. Since three implied slow relaxation timescales could be identified, we asked PCCA+ to cluster the trajectories into four metastable states (Figure 5). Moreover, PCCA+ provides the values of the MFPTs between each pair of metastable states (Figure 5), which were used to derive the rate constants corresponding to the entry of the solute into the membrane $k_{in}$, the exit out of it $k_{out}$, and the flip-flop transitions from one lipid layer to the other $k_{flip}$ (Figure 6A).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** PyEMMA implied timescales as a function of lag time for each Markov State Model. Timescales are sorted in descending order using the color lines blue, red, green, cyan, violet, yellow, and pink. Black lines indicate timescale equal to lag time.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** PyEMMA metastable states are represented by colored discs with a radius proportional to their population. The most populated states (labeled 3 and 4) always identify the water phase, while the smaller ones (labeled 1 and 2) represent the two lipid leaflets. Transitions are described with arrows whose thickness and label are associated to the rate constant and MFPT, respectively.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Designation of the kinetic rate constants in membrane planar (A) and spherical (B) geometry. In the liposome model, wo and lo stands for the outer water and lipid compartment, respectively, and similarly for the two inner compartments wi and li.</figcaption></figure>

## Kinetic Model Master Equations

Similarly to what was done by Dickson et al. [18], the master equations describing the passive diffusion of a solute from an outer water compartment into the water interior of a spherical liposome were solved by considering two lipid compartments (indicated as outer and inner lipid compartments) as schematically depicted in Figure 6B. The differential equations describing the time dependence of the solute in the four compartments read:

$$ \frac{dS_I(t)}{dt} = \sum_J [-k_{IJ}S_I(t) + k_{JI}S_J(t)] \qquad (5) \\ S_{wo}(0) = S_{tot}, \quad S_{lo}(0) = S_{li}(0) = S_{wi}(0) = 0 $$

where I, J = wo, lo, li, wi are the indices that refer to the outer water, outer lipid, inner lipid, and inner water compartments, respectively, and the summation runs over the nearest compartments. At time t = 0, the solute is entirely in the outer water compartment (Equation (5)).

Knowing that the rate constant k of each kinetic step depends on the permeability coefficient P, the surface area A separating two compartments, and the volume V of the donor compartment, through the relation k = P(A/V), it should be noted that the apparent rate constants k_{wo,lo}, k_{lo,wo}, k_{lo,li}, k_{li,lo}, k_{li,wi}, and k_{wi,li} in the membrane spherical geometry are related to the rate constants k_{in}, k_{out}, and k_{flip} previously estimated by the Markov State Models in the membrane planar geometry as

$$ \begin{aligned} k_{wo,lo} &= \frac{V_W}{V_{wo}} \frac{A_l}{A_L} k_{in}, & k_{wi,li} &= \frac{V_W}{V_{wi}} \frac{A_l}{A_L} k_{in} \\ k_{lo,wo} &= \frac{V_L}{V_l} \frac{A_l}{A_L} k_{out}, & k_{li,wi} &= \frac{V_L}{V_l} \frac{A_l}{A_L} k_{out} \\ k_{lo,li} &= \frac{V_L}{V_l} \frac{A_l}{A_L} k_{flip}, & k_{li,lo} &= \frac{V_L}{V_l} \frac{A_l}{A_L} k_{flip} \end{aligned} \qquad (6) $$

where A_L, V_L, and V_W are the lipid bilayer surface area, its interior volume, and the bulk water volume in the planar membrane systems simulated by MD, respectively. In the spherical liposome model, the surface area is considered equal on both outer and inner layers: A_l ≡ A_{li} = A_{lo}. Similarly, we assumed that the volumes of the two compartments in the lipid bilayer are equal: V_l ≡ V_{li} = V_{lo}. Quantitatively, we used a liposome radius r = 100 nm and a ratio of 1000:1:1:10 for the volumes of the four compartments V_{wo}: V_{lo}: V_{li}: V_{wi}, conforming to the experimental conditions of the liposome experiments reported in [10]. For the MD areas and volumes, we used the POPC area per lipid A_L = 65.6 Å² and the volume per lipid V_L = 1.2 nm³ [18]. The bulk
water volume $V_W$ was calculated from the volume per molecule 30.5 Å³ and the number of water per lipid (89 in our simulated systems).

Equations (5) are a system of four ordinary differential equations (three of which are independent) that we solved numerically using SciPy routines [38]. We also solved the system analytically, using the numerical values of the eigenvalues of the matrix associated with the system of differential equations (Appendix A).

## Permeability Coefficient from Kinetic Models

Following reference [39] and subsequent papers [10,40], the permeation coefficient is extracted from experimental data by fitting the luminescence curve with a biexponential function as a function of time:

$$I(t) = I(\infty) + [I(\infty)_I - I(0)_I]e^{-k_I t} + [I(\infty)_{II} - I(0)_{II}]e^{-k_{II} t} \quad (7)$$

where $I(\infty) = I(\infty)_I + I(\infty)_{II}$ is the total maximal luminescence, $I(\infty)_I$ and $I(\infty)_{II}$ are the maximal luminescence of the fast and slow phase respectively, $I(0)_I$ and $I(0)_{II}$ are the initial luminescence for the two phases, and $k_I$ and $k_{II}$ are the corresponding rate constants.

For luminescence experiments with liposomes, the permeation coefficient has been related to the rate constant $k_I$ of the fast exponential phase and the liposome radius as:

$$P_{lipo} = \frac{V_{wi}}{A_I} k_I = \frac{r}{3} k_I \quad (8)$$

From the point of view of the four-step model describing the process, where a peptide undergoes adsorption, flip-flop, and desorption through the liposome membrane bilayer, the relevant kinetics are expressed as three exponential functions, as in Equation (A7), rather than the two fitting functions for the luminescence intensity of Equation (7). We can, however, still use the rate constant associated with the fastest exponential phase in the liposome model and substitute its value in Equation (8) to obtain the membrane permeability coefficient, as it has been done in [18].

Errors on the permeability coefficient are derived using the Bayesian estimator in PyEMMA for the MSMs, which generates samples of transition matrices rather than one single transition matrix, as in the case of the maximum likelihood estimator. With this method, we obtained a sample of 100 transition matrices for each MSM, from which PyEMMA can compute samples for different quantities, such as the MFPT. The reported results were evaluated based on the mean value of the MFPT samples and their standard deviation, using propagation of uncertainties to compute the errors on the rate constants.

## 4. Conclusions

The ISDM approach is probably the most popular method used to estimate the membrane permeability of chemical compounds. However, Dickson et al. recently showed that an alternative method, based on Markov State Models of the permeation process, yielded permeability coefficients in better agreement with experimental data for seven small molecules [18]. In the present study, we applied and compared the two computational methods for a larger compound, the cyclic decapeptide CDP5. As opposed to the tendency reported by Dickson et al., we found that the ISDM-based and MSM-based method slightly underestimates and largely overestimates the peptide membrane permeability, respectively.

A probable explanation of the peptide permeability overestimation by our MSM-based calculations lies in the too-short duration of its unconstrained MD simulations required to build the Markov State Models. Indeed, due to its size and polarity, the CDP5 decapeptide cannot pass through the lipid headgroups as easily and rapidly as the small compounds simulated by Dickson et al. [18]. Thus, the peptide unconstrained simulations starting from maximum energy positions on the lipid headgroup surface in water require longer times to go down the hill and move towards the energy minimum inside the membrane. As a consequence, the peptide population $\pi$ and the free energy $\Delta G = -k_B T \log \pi$ at these maximum energy positions are larger and lower than expected,
respectively. Accordingly, the MSM-based peptide permeability coefficient is excessively large compared to experiments and ISDM calculations.

To remedy this default, it would probably be necessary to extend the duration of the unconstrained MD simulations and/or multiply their starting points. However, for the building of one Markov State Model, we have already performed a total of 15 µs MD simulations, whereas the ISDM approach has needed 10 µs US. In addition, the building of accurate MSMs requires running unconstrained MD simulations from maximum free energy states, which, very often, have to be predetermined by US calculations. Thus, predictions of the membrane permeability of peptides with the MSM-based method would require much more computer resources than the ISDM approach, which was able to yield reasonable estimations.

**Author Contributions:** Conceptualization, L.M. and T.H.-D.; Methodology, L.M. and T.H.-D.; Software, L.M. and T.H.-D.; Validation, L.M. and T.H.-D.; Formal analysis, L.M. and T.H.-D.; Investigation, L.M. and T.H.-D.; Resources, T.H.-D.; Data curation, L.M. and T.H.-D.; Writing—original draft, L.M. and T.H.-D.; Writing—review & editing, L.M. and T.H.-D.; Visualization, L.M. and T.H.-D.; Supervision, T.H.-D.; Project administration, T.H.-D.; Funding acquisition, T.H.-D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by HPC resources from GENCI-IDRIS (grants AD010713502 and A0130713807).

**Data Availability Statement:** All molecular dynamics trajectories supporting results reported in this study are available on request from the corresponding author. The raw data are not publicly available due to their size which prevent them from being uploaded to a publicly available repository.

**Conflicts of Interest:** The authors declare no conflict of interest regarding the publication of this article.

## Appendix A. Analytical Resolution of the Kinetic Model Equations

Equation (5) can be rewritten in matrix form as

$$ 
\begin{aligned}
\frac{d\vec{S}(t)}{dt} &= D \cdot \vec{S}(t) + \vec{B} \\
\vec{S}(0) &= \vec{S}_0
\end{aligned}
\qquad (A1)
$$

where $\vec{S}$ (and similarly $\vec{S}_0$) is defined as

$$ \vec{S} = \begin{pmatrix} S_{wo} \\ S_{-} \\ S_{wi} \end{pmatrix}, \quad S_{-} = \frac{1}{2}(S_{i} - S_{o}), \quad S_{i} = S_{wi} + S_{li} \qquad (A2) $$

The matrix $D$, the coefficients vector $\vec{B}$, and the initial conditions $\vec{S}_0$ read:

$$ D = \begin{pmatrix} -K_o & k_o & 0 \\ -k_f & -K_f & k'_f \\ 0 & k_i & -K_i \end{pmatrix} \qquad (A3) $$

$$ \vec{B} = \frac{S_{tot}}{2} \begin{pmatrix} k_o \\ k_{-} \\ k_i \end{pmatrix}, \quad \vec{S}_0 = \frac{S_{tot}}{2} \begin{pmatrix} 1 \\ -\frac{1}{2} \\ 0 \end{pmatrix} \qquad (A4) $$

where $S_{tot}$ indicates the constant total amount of solute in the system that is initially contained entirely in the outer aqueous compartment, as aforementioned. We have simplified notations by defining:

$$ 
\begin{aligned}
k_o &\equiv k_{lo,wo}, & K_o &\equiv k_o + k_{wo,lo} \\
k_i &\equiv k_{li,wi}, & K_i &\equiv k_i + k_{wi,li} \\
k_f &\equiv k_{li,lo}, & k'_f &\equiv k_{lo,li} \\
K_f &\equiv k_f + k'_f, & k_- &\equiv k'_f - k_f
\end{aligned}
\qquad (A5)
$$
## Appendix A.1. Approximated Solution for the Inner Acqueous Compartment

Knowing that the volume of the outer aqueous compartment is much larger than the volumes of the other compartments, and for symmetry reasons, we infer that

$$K_0 \approx k_0, \quad k_i \approx k_0, \quad K_f \approx 2k'_f \approx 2k_f \qquad (A6)$$

The solution to the differential equations yields the time evolution for the solute in the aqueous inner compartment, which is given by the sum of three exponential functions plus constants, in terms of the rate constants and of the eigenvalues $\lambda_1, \lambda_2, \lambda_3$ of the matrix $D$, with $\lambda_n = -\kappa_n$ and $\kappa_n \ge 0$, for $n = 1, 2, 3$:

$$S_{wi}(t) = \sum_{n=1}^{3} \left[ \frac{\beta_n}{\kappa_n} + \left( v_{0n} - \frac{\beta_n}{\kappa_n} \right) e^{-\kappa_n t} \right] \qquad (A7)$$

where $v_{0n}$ are the initial conditions for the eigenvectors $v_n$, and $\beta_n$ are the components of $\vec{B}$ in the eigenvector basis:

$$v_{0n} \equiv v_n(t = 0) = \frac{S_{tot}}{2} k_i \frac{K_i - \sum_{m \neq n} \kappa_m}{\prod_{m \neq n} (\kappa_m - \kappa_n)} \qquad (A8)$$

$$\beta_n = \frac{S_{tot}}{2} k_i \prod_{m \neq n} \frac{(K_i - \kappa_m)}{(\kappa_m - \kappa_n)}. \qquad (A9)$$

The initial conditions $v_{0n}$ and the coefficients $\beta_n$ are here derived in the approximations (A6), while the exact expressions are reported in Appendix A.2, together with the solutions for the solute in the remaining three compartments. In the limit $K_i \to k_i$, which implies $k_{wi,li} \ll k_{li,wi}$, the substrates are trivially constant as a function of time and equal to their initial values fixed in Equation (5):

$$S_{wo}(t) \equiv S_{tot}, \quad S_{lo}(t) \equiv S_{li}(t) \equiv S_{wi}(t) \equiv 0 \qquad (A10)$$

The condition $K_i \approx k_i$ would be satisfied if the inner aqueous compartment became much larger than the volume of the membrane, for instance, assuming that the rate constants are roughly of the same order of magnitude (as it will generally be the case for the current models). As we depart from this trivial limit, in other words, as the inner aqueous compartment volume becomes comparable to the membrane volume, we expect that the substrate would start decreasing in the outer aqueous compartment to partly fill the inner aqueous compartment once equilibrium is reached, as will be confirmed by our results.

## Appendix A.2. Full Solutions for the Kinetic Model Equations

The exact solutions to the set of Equation (5) describing the time evolution of the solute inside the liposome model illustrated in Figure 6B are sums of three exponential functions plus constants:

$$\begin{aligned}
S_{wo}(t) &= \sum_{n=1}^{3} \left[ \frac{k'_{f}}{k_{f}} - \frac{(K_{i} - \kappa_{n}) (k_{f} + k'_{f} - \kappa_{n})}{k_{i} k_{f}} \right] v_{n}(t) \\
S_{lo}(t) &= \sum_{n=1}^{3} \left[ \frac{k'_{f}}{k_{f}} - \frac{(K_{i} - \kappa_{n}) (k'_{f} - \kappa_{n})}{k_{i} k_{f}} \right] v_{n}(t) + \frac{S_{tot}}{2} \\
S_{li}(t) &= \sum_{n=1}^{3} \left[ \frac{k'_{f}}{k_{f}} - \frac{(K_{i} - \kappa_{n})}{k_{i}} \right] v_{n}(t) + \frac{S_{tot}}{2} \\
S_{wi}(t) &= \sum_{n=1}^{3} v_{n}(t)
\end{aligned} \qquad (A11)$$
where $v_n(t)$ are the eigenvectors solutions to the system of differential equations:

$$v_n(t) = \frac{\beta_n}{\kappa_n} + \left( v_{0n} - \frac{\beta_n}{\kappa_n} \right) e^{-\kappa_n t}, \quad n = 1, 2, 3 \qquad (A12)$$

If no specific approximation is applied for the membrane system, the initial conditions $v_{0n}$ for the eigenvectors $v_n$ yield:

$$v_{0n} = \frac{S_{tot}}{2} k_i \frac{K_i + k_- - \sum_{m \neq n} \kappa_m}{\prod_{m \neq n} (\kappa_m - \kappa_n)} \qquad (A13)$$

and the constant coefficients $\beta_n$ can be expressed as

$$\beta_j = \frac{S_{tot}}{2} k_i \frac{\prod_{m \neq n} (K_i - \kappa_m) + k_- (K_f + K_i - \sum_{m \neq n} \kappa_m) - k_f k_o + k'_f k_i}{\prod_{m \neq n} (\kappa_m - \kappa_n)} \qquad (A14)$$

## References

1. Ju, Y.; Guo, H.; Edman, M.; Hamm-Alvarez, S.F. Application of advances in endocytosis and membrane trafficking to drug delivery. *Adv. Drug Deliv. Rev.* **2020**, *157*, 118–141. [CrossRef] [PubMed]

2. Dean, M.; Moitra, K.; Allikmets, R. The human ATP-binding cassette (ABC) transporter superfamily. *Hum. Mutat.* **2022**, *43*, 1162–1182. [CrossRef] [PubMed]

3. Pizzagalli, M.D.; Bensimon, A.; Superti-Furga, G. A guide to plasma membrane solute carrier proteins. *FEBS J.* **2021**, *288*, 2784–2835. [CrossRef] [PubMed]

4. Vergalli, J.; Bodrenko, I.V.; Masi, M.; Moynié, L.; Acosta-Gutiérrez, S.; Naismith, J.H.; Davin-Regli, A.; Ceccarelli, M.; van den Berg, B.; Winterhalter, M.; et al. Porins and small-molecule translocation across the outer membrane of Gram-negative bacteria. *Nat. Rev. Microbiol.* **2020**, *18*, 164–176. [CrossRef]

5. Moynié, L.; Milenkovic, S.; Mislin, G.L.A.; Gasser, V.; Malloci, G.; Baco, E.; McCaughan, R.P.; Page, M.G.P.; Schalk, I.J.; Ceccarelli, M.; et al. The complex of ferric-enterobactin with its transporter from Pseudomonas aeruginosa suggests a two-site model. *Nat. Commun.* **2019**, *10*, 3673. [CrossRef]

6. Sugano, K.; Kansy, M.; Artursson, P.; Avdeef, A.; Bendels, S.; Di, L.; Ecker, G.F.; Faller, B.; Fischer, H.; Gerebtzoff, G.; et al. Coexistence of passive and carrier-mediated processes in drug transport. *Nat. Rev. Drug Discov.* **2010**, *9*, 597–614. [CrossRef]

7. Di, L.; Artursson, P.; Avdeef, A.; Ecker, G.F.; Faller, B.; Fischer, H.; Houston, J.B.; Kansy, M.; Kerns, E.H.; Krämer, S.D.; et al. Evidence-based approach to assess passive diffusion and carrier-mediated drug transport. *Drug Discov. Today* **2012**, *17*, 905–912. [CrossRef]

8. Lau, J.L.; Dunn, M.K. Therapeutic peptides: Historical perspectives, current development trends, and future directions. *Bioorg. Med. Chem.* **2018**, *26*, 2700–2707. [CrossRef]

9. Kansy, M.; Senner, F.; Gubernator, K. Physicochemical High Throughput Screening: Parallel Artificial Membrane Permeation Assay in the Description of Passive Absorption Processes. *J. Med. Chem.* **1998**, *41*, 1007–1010. [CrossRef]

10. Eyer, K.; Paech, F.; Schuler, F.; Kuhn, P.; Kissner, R.; Belli, S.; Dittrich, P.S.; Kramer, S.D. A liposomal fluorescence assay to study permeation kinetics of drug-like weak bases across the lipid bilayer. *J. Control. Release* **2014**, *173*, 102–109. [CrossRef]

11. Biedermann, F.; Ghale, G.; Hennig, A.; Nau, W.M. Fluorescent artificial receptor-based membrane assay (FARMA) for spatiotemporally resolved monitoring of biomembrane permeability. *Commun. Biol.* **2020**, *3*, 383. [CrossRef] [PubMed]

12. Di Pisa, M.; Chassaing, G.; Swiecicki, J.M. Translocation Mechanism(s) of Cell-Penetrating Peptides: Biophysical Studies Using Artificial Membrane Bilayers. *Biochemistry* **2015**, *54*, 194–207. [CrossRef]

13. Ulmschneider, J.P.; Ulmschneider, M.B. Molecular Dynamics Simulations Are Redefining Our View of Peptides Interacting with Biological Membranes. *Acc. Chem. Res.* **2018**, *51*, 1106–1116. [CrossRef] [PubMed]

14. Swift, R.V.; Amaro, R.E. Back to the Future: Can Physical Models of Passive Membrane Permeability Help Reduce Drug Candidate Attrition and Move Us Beyond QSPR? *Chem. Biol. Drug Des.* **2013**, *81*, 61–71. [CrossRef] [PubMed]

15. Kabelka, I.; Brozek, R.; Vacha, R. Selecting Collective Variables and Free-Energy Methods for Peptide Translocation across Membranes. *J. Chem. Inf. Model.* **2021**, *61*, 819–830. [CrossRef]

16. Marrink, S.J.; Berendsen, H.J.C. Simulation of water transport through a lipid membrane. *J. Phys. Chem.* **1994**, *98*, 4155–4168. [CrossRef]

17. Bennion, B.J.; Be, N.A.; McNerney, M.W.; Lao, V.; Carlson, E.M.; Valdez, C.A.; Malfatti, M.A.; Enright, H.A.; Nguyen, T.H.; Lightstone, F.C.; et al. Predicting a Drug’s Membrane Permeability: A Computational Model Validated with in Vitro Permeability Assay Data. *J. Phys. Chem. B* **2017**, *121*, 5228–5237. [CrossRef]

18. Dickson, C.J.; Hornak, V.; Pearlstein, R.A.; Duca, J.S. Structure–Kinetic Relationships of Passive Membrane Permeation from Multiscale Modeling. *J. Am. Chem. Soc.* **2017**, *139*, 442–452. [CrossRef]
19. Faulkner, C.; de Leeuw, N.H. Predicting the Membrane Permeability of Fentanyl and Its Analogues by Molecular Dynamics Simulations. *J. Phys. Chem. B* **2021**, *125*, 8443–8449. [CrossRef]

20. Sugita, M.; Sugiyama, S.; Fujie, T.; Yoshikawa, Y.; Yanagisawa, K.; Ohue, M.; Akiyama, Y. Large-Scale Membrane Permeability Prediction of Cyclic Peptides Crossing a Lipid Bilayer Based on Enhanced Sampling Molecular Dynamics Simulations. *J. Chem. Inf. Model.* **2021**, *61*, 3681–3695. [CrossRef]

21. Prinz, J.H.; Wu, H.; Sarich, M.; Keller, B.; Senne, M.; Held, M.; Chodera, J.D.; Schutte, C.; Noé, F. Markov models of molecular kinetics: Generation and validation. *J. Chem. Phys.* **2011**, *134*, 174105. [CrossRef]

22. Scherer, M.K.; Trendelkamp-Schroer, B.; Paul, F.; Pérez-Hernández, G.; Hoffmann, M.; Plattner, N.; Wehmeyer, C.; Prinz, J.H.; Noe, F. PyEMMA 2: A Software Package for Estimation, Validation, and Analysis of Markov Models. *J. Chem. Theory Comput.* **2015**, *11*, 5525–5542. [CrossRef]

23. Kumar, S.; Rosenberg, J.M.; Bouzida, D.; Swendsen, R.H.; Kollman, P.A. The weighted histogram analysis method for free-energy calculations on biomolecules. I. The method. *J. Comput. Chem.* **1992**, *13*, 1011–1021. [CrossRef]

24. Hub, J.S.; de Groot, B.L.; van der Spoel, D. g_wham-A Free Weighted Histogram Analysis Implementation Including Robust Error and Autocorrelation Estimates. *J. Chem. Theory Comput.* **2010**, *6*, 3713–3720. [CrossRef]

25. Wang, S.; König, G.; Roth, H.J.; Fouché, M.; Rodde, S.; Riniker, S. Effect of Flexibility, Lipophilicity, and the Location of Polar Residues on the Passive Membrane Permeability of a Series of Cyclic Decapeptides. *J. Med. Chem.* **2021**, *64*, 12761–12773. [CrossRef]

26. Fouché, M.; Schäfer, M.; Berghausen, J.; Desrayaud, S.; Blatter, M.; Piéchon, P.; Dix, I.; Martin Garcia, A.; Roth, H.J. Design and Development of a Cyclic Decapeptide Scaffold with Suitable Properties for Bioavailability and Oral Exposure. *ChemMedChem* **2016**, *11*, 1048–1059. [CrossRef] [PubMed]

27. Fouché, M.; Schäfer, M.; Blatter, M.; Berghausen, J.; Desrayaud, S.; Roth, H.J. Pharmacokinetic Studies around the Mono- and Difunctionalization of a Bioavailable Cyclic Decapeptide Scaffold. *ChemMedChem* **2016**, *11*, 1060–1068. [CrossRef] [PubMed]

28. Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.; Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; et al. CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. *J. Comput. Chem.* **2010**, *31*, 671–690. [CrossRef] [PubMed]

29. Huang, J.; MacKerell, A.D.J. CHARMM36 all-atom additive protein force field: Validation based on comparison to NMR data. *J. Comput. Chem.* **2013**, *34*, 2135–2145. [CrossRef]

30. O'Boyle, N.M.; Banck, M.; James, C.A.; Morley, C.; Vandermeersch, T.; Hutchison, G.R. Open Babel: An open chemical toolbox. *J. Cheminform.* **2011**, *3*, 33. [CrossRef]

31. Schrödinger, LLC. *The PyMOL Molecular Graphics System*, Version 1.8; Schrödinger, LLC.: Cambridge, MA, USA, 2015.

32. Essmann, U.; Perera, L.; Berkowitz, M.L.; Darden, T.; Lee, H.; Pedersen, L.G. A smooth particle mesh Ewald method. *J. Chem. Phys.* **1995**, *103*, 8577–8593. [CrossRef]

33. Hess, B.; Bekker, H.; Berendsen, H.J.C.; Fraaije, J.G.E.M. LINCS: A linear constraint solver for molecular simulations. *J. Comput. Chem.* **1997**, *18*, 1463–1472. [CrossRef]

34. Lee, J.; Cheng, X.; Swails, J.M.; Yeom, M.S.; Eastman, P.K.; Lemkul, J.A.; Wei, S.; Buckner, J.; Jeong, J.C.; Qi, Y.; et al. CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field. *J. Chem. Theory Comput.* **2016**, *12*, 405–413. [CrossRef]

35. Abraham, M.J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J.C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. *SoftwareX* **2015**, *1–2*, 19–25. [CrossRef]

36. Jorgensen, W.L.; Chandrasekhar, J.; Madura, J.D.; Impey, R.W.; Klein, M.L. Comparison of simple potential functions for simulating liquid water. *J. Chem. Phys.* **1983**, *79*, 926–935. [CrossRef]

37. Röblitz, S.; Weber, M. Fuzzy spectral clustering by PCCA+: Application to Markov state models and data classification. *Adv. Data Anal. Classif.* **2013**, *7*, 147–179. [CrossRef]

38. Virtanen, P.; Gommers, R.; Oliphant, T.E.; Haberland, M.; Reddy, T.; Cournapeau, D.; Burovski, E.; Peterson, P.; Weckesser, W.; Bright, J.; et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nat. Methods* **2020**, *17*, 261–272. [CrossRef] [PubMed]

39. Thomae, A.V.; Wunderli-Allenspach, H.; Krämer, S.D. Permeation of aromatic carboxylic acids across lipid bilayers: The pH-partition hypothesis revisited. *Biophys. J.* **2005**, *89*, 1802–1811. [CrossRef]

40. Thomae, A.V.; Koch, T.; Panse, C.; Wunderli-Allenspach, H.; Krämer, S.D. Comparing the lipid membrane affinity and permeation of drug-like acids: The intriguing effects of cholesterol and charged lipids. *Pharm. Res.* **2007**, *24*, 1457–1472. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.