# Probing Passive Permeation of Tetracycline: Are Simulations Ready for beyond-Rule-of-Five Drug Permeability Calculation?

Yajing Qi, Christophe Chipot, Yi Wang

► To cite this version:

Yajing Qi, Christophe Chipot, Yi Wang. Probing Passive Permeation of Tetracycline: Are Simulations Ready for beyond-Rule-of-Five Drug Permeability Calculation?. Journal of Physical Chemistry B, 2025, 129 (41), pp.10810-10823. 10.1021/acs.jpcb.5c05445. hal-05412281

HAL Id: hal-05412281

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

Submitted on 11 Dec 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
Probing Passive Permeation of
beyond-Rule-of-Five Drugs: Are Simulations Up
to the Task?

Yajing Qi ID,† Christophe Chipot ID,*,‡ and Yi Wang ID*,†

†Department of Physics, The Chinese University of Hong Kong, Shatin, Hong Kong SAR, China

‡Laboratoire International Associé Centre National de la Recherche Scientifique et University of
Illinois at UrbanaChampaign, Unité Mixte de Recherche n°7019, Université de Lorraine, 54506

Vandœuvre-lès-Nancy, France

¶Theoretical and Computational Biophysics Group, Beckman Institute, and Department of

Physics, University of Illinois at UrbanaChampaign, Urbana, IL 61801, USA

§Department of Biochemistry and Molecular Biology, University of Chicago, Chicago, IL 60637,

USA

E-mail: Christophe.Chipot@univ-lorraine.fr; yiwang@cuhk.edu.hk
Abstract

Passive permeation across lipid membranes is a key determinant of drug bioavailability and efficacy. Accurate computational estimation of drug permeability is essential for ratio-
nal drug design, yet remains challenging, particularly for ionizable and beyond-Rule-of-Five
(bRo5) compounds. In this study, we employ advanced molecular simulations and the inho-
mogeneous solubility-diffusion model to calculate the effective permeability and elucidate the
membrane permeation mechanism of the antibiotic tetracycline (TC), the six hydrogen-bond
donors of which violates one of Lipinski's Rule-of-Five. By integrating the pH-partitioning
and Boltzmann-weighted average potential schemes and accounting for both its neutral (TCN)
and zwitterionic (TCz) tautomers, we show that the dominant contribution to the effective per-
meability of TC arises from TCN, despite its low abundance. This result is attributed to the
relatively small permeation barrier of TCN, explaining why the antibiotic exhibits moderate
effective permeability even though it predominantly exists in the zwitterionic form at neutral
pH. A further systematic investigation reveals that membrane patch size significantly impacts
permeability estimates for TC, in contrast to the relative insensitivity observed for three other
permeants. This unique sensitivity can be attributed to the hydrogen-bond network formed
between TC and its lipid environment, with the smallest 32-POPC patch artificially raising
the drug molecule's permeation barrier and the largest 256-POPC patch introducing signifi-
cant hysteresis. Overall, our results suggest that while probing the passive permeation of bRo5
drugs by molecular simulations appears increasingly feasible, the protonation and tautomeric
states of the permeants, the estimation of their microscopic acid dissociation constants, as well
as the potential impact of membrane patch-size effects need to be fully considered in perme-
ability predictions for these complex molecules.
# Introduction

The passive permeation of drug molecules across lipid membranes critically influences their bioavailability and therapeutic efficacy.¹ Gaining molecular-level insight into how small molecules traverse lipid bilayers and quantifying the associated membrane permeability is key to the rational design of drugs with desirable ADMET (adsorption, distribution, metabolism, excretion, and toxicology) properties.² Experimentally, such permeability is often estimated from the flux across a monolayer of cells (Caco-2) or through the parallel artificial membrane permeability assay (PAMPA).³⁻⁵ Computationally, molecular dynamics (MD) simulations are frequently employed to investigate the thermodynamic and kinetic underpinnings of passive permeation.⁶⁻⁹ While alternative approaches have yielded notable successes,¹⁰,¹¹ the majority of these simulations have been based on the “inhomogeneous solubility-diffusion” (ISD) model,¹²,¹³ which combines the potential of mean force (PMF) underlying the translocation of a permeant molecule across a lipid bilayer with its position-dependent diffusivity profile to yield the membrane permeability (P):

$$ \frac{1}{P} = \int_{z_1}^{z_2} dz \frac{e^{+\beta w(z)}}{D(z)}, \qquad (1) $$

where the membrane is assumed to span over $[z_1, z_2]$ with its normal placed along the z-axis; The PMF $w(z)$ measures the free energy profile as a function of the permeant molecule’s z-coordinate, while $D(z)$ stands for its position-dependent diffusion coefficient; $\beta = (k_B T)^{-1}$, with $k_B$ representing the Boltzmann constant and $T$ the temperature. Based on Eq. 1, permeability calculations from molecular simulations have been performed on a variety of drug or drug-like small molecules.¹⁰,¹¹,¹⁴⁻¹⁹ Deviations of calculated permeabilities from the corresponding experimental references have been found to be as low as 0.12 log unit,¹⁹ although deviations greater than one log unit are frequently reported.⁹,¹¹,¹⁷,²⁰ A major contribution to errors in the computed permeabilities is the PMF calculation—since $w(z)$ appears in the exponent of Eq. 1, even small errors in the PMF can be amplified, resulting in significant uncertainty in the calculated permeability.

For ionizable molecules, an additional complexity arises in their permeability calculation:
With the exception of constant-pH treatment,²¹ molecular simulations are conducted for a single ionization state of the permeant. For example, an acid molecule with one ionizable group (HA ⇌ H⁺ + A⁻) has two ionization states. A molecular simulation containing either the negatively charged species A⁻ or the neutral species HA can only produce the PMF and diffusivity profiles for that species. The respective permeability obtained via Eq. 1, termed specific permeability,¹⁸ reflects the intrinsic membrane-crossing capacity of the permeant in the given ionization state, and, by definition, is pH-independent. The combined permeability of all ionization states of a given molecule, termed the effective permeability, is pH-dependent and needs to be estimated from the PMF and diffusivity calculations of its individual ionization states. Recently, Harris et al.¹⁸ thoroughly evaluated two such protocols, namely, the pH-partitioning and the Boltzmann-weighted average potential (BWAP). In the former scheme, following the pH-partition hypothesis, the permeant molecules in bulk solution are assumed to be an equilibrium mixture of each ionization state, with the fraction of a given state determined by the Henderson-Hasselbalch equation. P_s, the specific permeability for each ionization state s, is calculated using Eq. 1 from the corresponding PMF and diffusivity, w_s(z) and D_s(z), respectively, and then combined to yield the effective permeability:¹⁸

$$P_{\text{eff}} = \sum_{s} f_{s} P_{s}, \qquad (2)$$

where f_s represents the fraction of the ionization state s at the given pH. The pH-partitioning protocol assumes that the protonation/deprotonation transitions of the permeant are much slower than its diffusion along z, i.e., the molecule permeates the membrane before it has time to switch its ionization state. In contrast, permeability calculation via BWAP, which underlies the constant-pH treatment,²¹ assumes that the protonation/deprotonation transitions are much faster than diffusion. As a result, an equilibrium between each ionization state is always established along z, and the permeant diffuses on a BWAP, w_m(z), given by:¹⁸

$$e^{-\beta w_m(z)} = \sum_{s} f_s e^{-\beta w_s(z)}. \qquad (3)$$
This BWAP, along with an effective diffusion coefficient (Eq. 17 of Harris, et al. ¹⁸), is then used to compute $P_{\text{eff}}$ of the ionizable permeant via Eq. 1. In general, whether the pH-partitioning or BWAP protocol should be employed to calculate $P_{\text{eff}}$ depends on how fast the permeant crosses the membrane relative to its protonation/deprotonation. For drug molecules with a high effective permeability, such as nicotine and varenicline at neutral pH, the pH-partitioning scheme has been found to better characterize the physics of their permeation processes. ¹⁸

The importance of considering protonation/deprotonation kinetics has also been discussed recently by Sezer and Oruç²² in their analysis of liposomal fluorescence assays designed by Eyer et al.²³ to measure the permeability of drug molecules that are either weak acids or weak bases. Compared with PAMPA and Caco-2 assays, these liposome-based assays have negligible unstirred water layer (UWL) effect²⁴ and no paracellular contributions, making the measured permeability easier to interpret.¹⁸ In addition, they utilize single-component bilayers (e.g., POPC), allowing the measured *apparent* permeability ($P_{\text{app}}$) to be directly compared with the computed $P_{\text{eff}}$ from molecular simulations conducted with the same lipid composition.

Based on liposomal fluorescence assays, Kramers et al.²⁵ and Hermann et al.²⁶ investigated the permeability of a series of drug molecules, including, notably, some beyond-Rule-of-Five (bRo5) ones. Formulated by Lipinski and colleagues in 1997,²⁷ the Rule-of-Five (Ro5) predicts that poor absorption or permeation is more likely when a molecule has one or more violations of the following: no more than five hydrogen bond donors, no more than ten hydrogen bond acceptors, a molecular weight less than 500 and a calculated logarithm of partition coefficient (log P) no greater than five. The Ro5 was set such that among the thousands of drugs in the USAN library investigated by Lipinski et al., there were about 10% outliers for each of the above rules.²⁷ The percentage of breaking the “no more than five hydrogen bond donors” rule has been found to be the lowest,²⁵ and orally-dosed drugs approved over the past several decades each have ~2 hydrogen bond donors on average.²⁸ While compounds violating one or more Ro5 are generally expected to have poor absorption or permeation, exceptions do exist—notably the widely used antibiotic tetracycline (TC). TC and its analogs can treat a broad range of bacterial infections,²⁹ some of which have also been
shown to inhibit the growth of drug-resistant *M. tuberculosis*.³⁰ With six hydrogen-bond donors, TC has one Ro5 violation. The drug molecule has four fused rings and three ionizable groups (Fig. 1A, B). Using the aforementioned liposomal assays, Kramer et al.²⁵ demonstrated that TC has a moderate permeability at pH = 6 ($P_{app} = 10^{-5.86}$ cm/s). At this pH, TC primarily exists as a zwitterion,³¹ i.e., the molecule carries a zero net charge, but contains both negatively charged and positively charged groups. While the measured $P_{app}$ is consistent with previous results of a PAMPA assay, as well as the fast adsorption and good bioavailability of the drug *in vivo*,³²,³³ the charged groups carried by the zwitterionic drug have been suggested to significantly reduce its membrane permeability at neutral pH.³⁴ Thus, the origin of the moderate membrane permeability to TC remains unclear.

<figure><img src="image_2.png" /><figcaption>(C) Phenolic diketone</figcaption></figure>

<caption>(B) Dimethylammonium</caption>

<caption>(A) Tricarbonyl</caption>

<figure><img src="image_3.png" /><figcaption>C. Methanol</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>D. Urea</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>E. 2',3'-Dideoxyadenosine (DDA)</figcaption></figure>

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

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

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

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

**Figure 1:** Structural representations of the permeants investigated in this work. The top part of each panel (A–E) shows their chemical formula, while the bottom part shows their three-dimensional structures rendered with ±5 $k_B T/e$ electrostatic potential isosurfaces based on atomic partial charges.
In this work, we set out to explore the membrane permeation mechanism of TC through computing its effective permeability. At pH = 6, two important tautomers of the zero-net-charge TC are the uncharged neutral form (TC<sub>N</sub>) and the zwitterionic form (TC<sub>Z</sub>), the latter of which is obtained from the former through transferring a proton from the tricarbonyl group A to the dimethylammonium group B (Fig. 1 A, B). Through determining the PMFs and diffusivity profiles for both TC<sub>Z</sub> and TC<sub>N</sub>, we show that unlike previously hypothesized,<sup>35</sup> the positively charged dimethylammonium and the negatively charged tricarbonyl of TC<sub>Z</sub> do not sufficiently cancel each other's charge to produce a significant membrane permeability. Instead, the effective permeability of TC in a wide range of pH values is predominantly attributed to the TC<sub>N</sub> tautomer, despite its much smaller population than TC<sub>Z</sub>. Only by considering both zwitterionic and neutral tautomers in the framework of pH-partitioning or BWAP, a properly weighted P<sub>eff</sub> can be obtained, which falls within approximately one log unit of P<sub>app</sub> measured experimentally.

With its six hydrogen-bond donors, significant membrane deformation is recorded during the permeation of TC. This observation further prompts us to conduct a systematic evaluation of bilayer patch size on the membrane permeability calculation of TC as well as three other permeants varying in their size and polarity, namely, methanol, urea and 2',3'-dideoxyadenosine (DDA). The influence of the patch size on membrane properties—including lipid diffusion, bilayer structure, and elastic moduli—has been well documented.<sup>36-38</sup> Smaller membrane patches, often necessitated by computational constraints, can artificially suppress long-wavelength fluctuations, alter lateral pressure profiles, decrease lipid and peptide diffusion, thereby, potentially biasing the permeability estimates. Despite these previous studies, the sensitivity of permeability calculation to bilayer patch size remains underexplored, particularly for drug or drug-like molecules. Through comparative simulations in POPC bilayers with 32, 64, 128 or 256 lipids, we quantify how patch size modulates the thermodynamic barriers associated with membrane permeation of the aforementioned four molecules. Notably, we find that only TC demonstrates a strong patch size dependence in its computed permeability, where the smallest 32-POPC bilayer presents artificially large barriers against its permeation, while the biggest patch with 256 POPC shows sampling hysteresis that
compromises the quality of the free-energy calculation. We attribute this sensitivity to the strong
hydrogen-bonding capacity of TC, a hallmark of its Ro5 violation.

Collectively, our results suggest that while permeability calculations for bRo5 drugs based on
molecular simulations and the ISD model are increasingly feasible, care must be taken to ensure
that (1) for ionizable molecules that can exist in both zwitterionic and neutral forms, both tautomers
are simulated and their results combined with other ionization states using the pH-partitioning or
BWAP protocols, and (2) a bilayer patch size that is neither too small (to avoid artificially suppress-
ing membrane deformation) nor too big (to minimize hysteresis) is adopted in these simulations.
Below, we describe our computational methods and results in detail, followed by a discussion on
challenges facing bRo5 drug permeability calculation. We note that apart from errors in free-energy
and diffusivity calculations, uncertainty in estimating microscopic p*K<sub>a</sub> values can contribute sig-
nificantly to deviations between computed P<sub>eff</sub> and experimentally measured P<sub>app</sub>. In particular,
this issue is exacerbated if incorrectly assigned macroscopic p*K<sub>a</sub> values are used in place of the
appropriate microscopic ones.

Methods

All MD simulations reported in this work were performed in 1-palmitoyl-2-oleoyl-sn-glycero-3-
phosphocholine (POPC) bilayers. Unless otherwise specified, a POPC bilayer with 128 lipids was
used by default. Three other bilayer patch sizes, namely, 32, 64, and 256 POPC, were subse-
quently used to examine the effect of patch size on permeability calculation. All bilayer systems
were built using the CHARMM GUI program.³⁹ For bilayers with 32, 64, 128 and 256 POPC, the
total number of atoms was approximately 12,000, 21,000, 39,000, and 81,000, respectively, while
the system sizes were approximately 34 × 34 × 105 Å³, 48 × 48 × 95 Å³, 66 × 66 × 90 Å³, and
95 × 95 × 90 Å³, respectively. All MD simulations were carried out using NAMD 3.0,⁴⁰,⁴¹ with
the CHARMM36 force field (FF) for lipids and TIP3P for water.⁴²,⁴³ The force field parameters
for methanol, urea, and 2',3'-dideoxyadenosine were obtained from the CGenFF web server.⁴⁴,⁴⁵
Parameters of TC were taken from Aleksandrov and Simonson,⁴⁶,⁴⁷ in which parameters compatible with the CHARMM FF were developed for tetracycline and 11 of its analogs. The system temperature was maintained at 310 K using a Langevin thermostat with a damping coefficient of 1 ps⁻¹,⁴⁸ with the pressure regulated at 1.01325 bar via the Langevin piston method.⁴⁹ Long-range electrostatic interactions were evaluated employing the particle-mesh Ewald algorithm with a grid spacing of 1 Å,⁵⁰ whereas short-range van der Waals (vdW) and electrostatic interactions were truncated at a smoothed 12 Å cutoff. The vdW potential was treated with a force-switching protocol. Multiple time-stepping was used to integrate the equations of motion, with time-steps of 2 fs for short-range and 4 fs for long-range interactions.⁵¹ Water molecules were constrained to their equilibrium geometry using the SETTLE algorithm,⁵² and covalent bonds involving hydrogen atoms were constrained with the RATTLE algorithm.⁵³

To determine the free-energy profile associated with the membrane permeation of a given molecule, we employed the Well-Tempered Metadynamics- extended-system Adaptive Biasing Force (WTM-eABF) algorithm⁵⁴,⁵⁵ with the collective variable (CV) defined as the z-projection of the center of mass of the permeant. The center of mass of lipid phosphorus atoms was placed at z = 0. The overall range of sampling, set to −43 < z < 43 Å (−38 < z < 38 Å for systems with 256 POPC), was divided into three windows. Each window spanned 28 to 33 Å, with a 5 Å overlap between adjacent windows. During the WTM-eABF calculation, sampling data were accumulated and recorded using a bin width of 0.2 Å for each window. To minimize spurious nonequilibrium phenomena, the biasing potential was withheld until 20,000 sampling steps had been collected in a given bin.⁵⁶ For the WTM component, the bias temperature was set to 3000 K, with a Gaussian hill width of 5.0 Å and a hill height of 0.1 kcal/mol. For the eABF component, the standard deviation between the colvar and the fictitious particle was 0.05 Å, with an oscillation period of 100 ns.

After obtaining a well-converged PMF, we evaluated the diffusivity of a permeant by applying its inverse PMF in the *z* direction. As a result of the applied inverse PMF, the permeant was able to diffuse freely over an approximately flat free-energy surface along the membrane normal. This was achieved by initializing an ABF calculation with the final, converged PMF, disabling further
updates of the biasing force. Since diffusion in the (x, y)-plane was inherently unrestrained, introducing the reverse biasing force in the z-direction allowed for free diffusion in the entire simulation box. The position-dependent fractional diffusivity, as well as the classical diffusivity, were estimated using a variant of the Bayesian inference scheme developed for classical diffusion.⁵⁷ These calculations were performed with the DiffusionFusion program,⁵⁸,⁵⁹ employing a grid spacing of 1 Å. Finally, analysis of permeant orientation, hydrogen bonds, and the number of water molecules accompanying the permeants in the region $-15 \AA \le z \le 15 \AA$ was performed using VMD.⁶⁰ The average values and standard deviations for each property were computed for every 1-Å bin along the z-axis. Hydrogen bonds were determined using a cutoff distance of 3.6 Å and an angle of 30°. The number of water within 2.4 Å of the permeants was also computed.

# Results

## The specific membrane permeabilities to TC<sub>N</sub> and TC<sub>Z</sub> differ by six orders of magnitude

We begin our investigation of TC membrane permeation by computing the free-energy profiles of two zero-net-charge TC tautomers (Fig. 1) crossing a bilayer with 128 POPC lipids, using the WTM-eABF approach. Given that all bilayers in our simulations contain the same number and type of lipids in its two monolayers, the resulting PMFs are symmetrized (see Methods), with the differences between the original, un-symmetrized PMFs and the final, symmetrized ones shown as error bars in Fig. 2. All simulations are continued until the corresponding error bars are reduced to approximately 1 kcal/mol or less. The sampling quality is further evaluated based on the degree of convergence of the computed PMFs and their gradients: As shown in Fig. S1, the PMFs obtained with increasing simulation time offer a direct assessment of the evolution of the free-energy profiles, where sampling convergence is reflected in the small ($< 0.5$ kcal/mol) differences between the PMF profiles computed using the last ~0.5-μs simulations. With the final PMF gradient as reference, we also compute the root mean square deviation (RMSD) of the PMF gradients obtained with increasing simulation time, the evolution of which (Fig. S1) further confirms
the convergence of our free-energy calculations. Next, the diffusivity profiles of the two TC tau-
tomers are determined from additional sets of 4-μs and 2-μs simulations with their inverse PMFs
applied as biasing potentials, respectively. The permeants experience an approximately flat free
energy surface in these simulations, enabling them to undergo free diffusion along the membrane
normal. The position-dependent classical diffusivity profiles computed from the simulation trajec-
tories are illustrated in Fig. 2B, while the fractional diffusivity profiles are shown in Fig. S2. As
with the PMFs, the error bars in Fig. 2B and Fig. S2 represent the differences between the original,
un-symmetrized diffusivity profiles and the ones after symmetrization. It is noteworthy that the
diffusivity profiles obtained with the classical and fractional Smoluchowski equations are nearly
identical, revealing no hallmark of anomalous diffusion,⁶¹ which is suggestive that the issue of
timescale separation of the CV (projection along the z-axis) and the slow degrees of freedom of
nearby lipids—or the incompleteness of the CV subspace describing the permeation event—may
be less acute for large permeants.¹⁷,⁵⁸ Finally, using the ISD model as given in Eq. 1, the specific
permeability of the neutral TC<sub>N</sub> (P<sub>N</sub>) and the zwitterionic TC<sub>Z</sub> (P<sub>Z</sub>) are determined from their
respective w(z) and D(z) profiles: P<sub>N</sub> = 2.22 × 10<sup>-2</sup> cm/s and P<sub>Z</sub> = 1.63 × 10<sup>-8</sup> cm/s. From
hereon, we frequently refer to the 10-based log values of the permeabilities (log P<sub>N</sub> = -1.65 and
log P<sub>Z</sub> = -7.79) with their units (cm/s) omitted for simplicity.

Our results show that P<sub>N</sub> is approximately six orders of magnitude greater than P<sub>Z</sub>, indicating
that while both tautomers carry zero net charge, the neutral TC<sub>N</sub> is considerably more membrane
permeable than the zwitterionic TC<sub>Z</sub>. As their diffusivity profiles are nearly identical, the large dis-
crepancy in their specific permeabilities arises almost exclusively from their distinct permeation
free energy profiles. While TC<sub>N</sub> and TC<sub>Z</sub> both have “W-shaped” PMFs, as observed in the mem-
brane translocation PMFs of many other small molecules,<sup>6,15,20,62</sup> the permeation energy barriers
encountered by the two tautomers differ significantly: the zwitterionic TC<sub>Z</sub> faces a barrier of 13.6
kcal/mol at the bilayer center (z = 0), while the neutral TC<sub>N</sub> has a barrier of only 4.5 kcal/mol at
the same location. Since w(z) appears in the exponent of the integrand in Eq. 1, the 9.1 kcal/mol
difference in their energy barriers results in the six orders-of-magnitude difference in the specific
permeabilities of the two tautomers.

<figure><img src="image_10.png" /><figcaption>Figure 2: (A) Symmetrized PMFs $w(z)$ and (B) classical diffusivity profiles $D(z)$ of TC$_N$ and TC$_Z$. PMFs for TC$_N$ both with and without the free energy shift $\Delta G = -k_B T \ln(f_N/f_Z)$, are shown. (C) Complete ionization scheme of TC. Black arrows correspond to microscopic $pK_a$s and blue arrows indicate the transition between neighboring tautomers with the same net charge. One thermodynamic cycle involving $A^-B^+C^0$, $A^-B^0C^0$ and $A^0B^0C^0$ is highlighted.</figcaption></figure>

To understand the origin of their distinct permeation free energy profiles, we analyzed the orientations of TC$_N$ and TC$_Z$ using two vectors: $\vec{v}_1$ pointing from atom C2 near the tricarbonyl group to atom C11 on the aromatic ring D, and $\vec{v}_2$ pointing from atom C15 near the phenolic dike-tone group to atom C22 of the methyl group attached to ring C (Fig. 3 A). The position-dependent profiles of their angles with the z-axis reveal distinct orientations of the two tautomers (Fig. 3 B, C): For TC$_N$, $\vec{v}_1$ maintains an approximately 90° angle to the z-axis, while the angle between $\vec{v}_2$ and the z-axis switches from ~150° to ~30° as TC$_N$ crosses the bilayer center, indicating a flip-flop transition along $\vec{v}_2$. In contrast, for TC$_Z$, $\vec{v}_2$ maintains an approximately 90° angle with the membrane normal, while a flipflop transition is observed along the vector $\vec{v}_1$. Together with representative snapshots from their simulations shown in Fig. 4, these data indicate that while both
tautomers undergo a flipflop transition, they have distinct orientational preferences: TC<sub>N</sub> prefers an approximately “horizontal” orientation, placing its four fused rings perpendicular to the membrane normal and exposing its phenolic diketone and amide groups to water; TC<sub>Z</sub> adopts a predominantly “vertical” orientation and exposes its charged dimethylammonium and tricarbonyl to water. These behaviors can be explained by the distinct polarity profiles of the two tautomers: The deprotonation of the dimethylammonium significantly reduces the polarity of this moiety, rendering the TC<sub>N</sub> tautomer nearly amphiphilic—four of its six hydrogen bond donors are placed on the phenolic diketone side and of the remaining two hydrogen bond donors, the hydroxyl from the tricarbonyl group readily forms an intramolecular hydrogen bond. Consistent with its polarity profile, the TC<sub>N</sub> tautomer experiences a deep energy well (−4.9 kcal/mol) at the lipid-water interface, where it can expose its highly polar phenolic diketone side to water and bury the opposite, less polar side in the lipid tails. The TC<sub>Z</sub> tautomer, in contrast, strongly favors exposing its positively charged dimethylammonium and negatively charged tricarbonyl to water. Since both of these charged components are near ring A (Fig. 3 A), the “vertical” orientation that maximizes their exposure to water can only be achieved at the expense of burying the polar phenolic diketone group in lipid tails, which results in only a shallow energy well (−2.1 kcal/mol) for TC<sub>Z</sub> at the lipid-water interface. The zwitterionic components in TC<sub>Z</sub> also appear to elicit distinct hydrogen bond interactions and solvation shells from that of TC<sub>N</sub>, with the former solvated by ~6–8 water molecules throughout its membrane permeation, compared to ~3–5 water molecules accompanying the latter tautomer (Fig. 3 D-H). Overall, the charges carried by its zwitterionic components and the resulting greater requirement for solvation underlies the significantly higher energy barrier against the permeation of TC<sub>Z</sub> compared to the neutral TC<sub>N</sub>.

**The effective permeability of TC at pH = 6 is dominated by TC<sub>N</sub>**

To compute the effective permeability of TC at the weakly acidic condition (pH = 6) used in liposomal fluorescence experiments,²⁵,²⁶ we adopt the pH partitioning and BWAP protocols as given in Harris et al. ¹⁸ Apart from membrane permeation PMFs and diffusivity profiles, both protocols
<figure><img src="image_11.png" /></figure>

Figure 3: Orientation, the number of solvating water molecules and hydrogen bond (H-bond) interactions of TC during its membrane permeation. (A) Definition of $\vec{v}_1$ and $\vec{v}_2$. (B, C) The angle between $\vec{v}_1$ (B) or $\vec{v}_2$ (C) and the z-axis. (D) The number of water molecules within 2.4 Å of TC. (E-H) The number of total H-bonds (E), TC's intramolecular H-bonds (F), H-bonds between TC and POPC lipids (G) or water (H). Error bars represent the standard deviations.

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

Figure 4: Representative snapshots of TC<sub>N</sub> and TC<sub>Z</sub> during membrane permeation. Hydrogen bonds within TC as well as those between TC and POPC or water molecules are shown in black dashed lines.
require the ratios, or equivalently, the free-energy differences between various protonation mi-
crostates of the permeant molecule. As explained by Gunner et al.,⁶³ a protonation macrostate is
defined by the net charge of the molecule, while a microstate is defined by the specific protonation
and tautomeric states of all its ionizable groups: For a simple molecule with a single ionizable
group, each macrostate may contain only one microstate; for a complex molecule with multiple
ionizable groups, a macrostate may contain multiple microstates that represent different tautomers
of the molecule at the given net charge. With three ionizable groups, TC has four protonation
macrostates characterized by a net charge of +1, 0, -1 and -2, respectively. Its complete ioniza-
tion scheme showing all eight microstates is given in Fig. 2 C, where the superscripts indicate the
charge carried by the tricarbonyl (A), the dimethylammonium (B) and the phenolic diketone group
(C), respectively (Fig. 1 A). For instance, A⁰B⁰C⁰ represents TCₙ, and A⁻B⁺C⁰ corresponds to
its tautomer TCz. The remaining zero-net-charge tautomer A⁰B⁺C⁻ is also a zwitterion and will
be referred to as TCz⁻ from hereon. In theory, both the pH partitioning protocol given in Eq. 2 and
the BWAP treatment given in Eq. 3 should enumerate all microstates of a given molecule. In prac-
tice, calculating the membrane permeation free energy and diffusivity profiles for all microstates
is neither feasible nor necessary, since only those microstates with significant fraction and/or high
specific permeability make non-negligible contributions to Pₐff. At pH = 6, the fractions of the four
macrostates of TC can be readily obtained from its three macroscopic pKₐs, consistently found to
be around 3, 7, and 9, respectively.³¹,⁶⁴,⁶⁵ In this work, we adopt the values reported by Leeson
et al.,⁶⁴ namely, pKₐ₁ = 3.33, pKₐ₂ = 7.75 and pKₐ₃ = 9.61, which are very similar to those
reported by Stephens et al.³¹ and Qiang and Adams⁶⁵ (see Table S1). As shown in Table 1, with
a fraction of ~98 %, the zero-net-charge macrostate dominates at pH = 6, while the fractions of
the +1 and -2 charged macrostates are three and six orders of magnitude smaller, respectively.
The -1 charged macrostate is the second most populated state at this pH, which has a fraction of
~1.74%. Since our calculations have revealed that the ±1 charged components in TCz drastically
reduce its specific permeability, all charged macrostates of the antibiotic can be expected to have
specific permeabilities at least as low as the zwitterionic TCz. Therefore, in subsequent calcula-
tions, we assume that their specific permeabilities are capped by that of TC<sub>Z</sub> (or equivalently, their membrane permeation PMFs lower-bounded by that of TC<sub>Z</sub>). Similarly, the alternative zwitterionic tautomer TC<sub>Z</sub>' is assumed to have the same specific permeability as TC<sub>Z</sub>. As will be shown next, this assumption has negligible impact on the computed value of P<sub>eff</sub>, as long as the neutral TC<sub>N</sub> remains significantly more permeable.

<table><caption>**Table 1:** The macrostates of TC and their fractions at pH = 6. Instead of f<sub>2</sub>, fractions of the TC<sub>Z</sub>, TC<sub>N</sub>, and TC<sub>Z'</sub> tautomers are shown.</caption><thead><tr><th>Macrostate</th><th>s = 3</th><th>s = 2</th><th>s = 1</th><th>s = 0</th></tr></thead><tbody><tr><th>Net charge</th><td>+1</td><td>0</td><td>-1</td><td>-2</td></tr><tr><th>f<sub>s</sub></th><td>2.10 &times; 10<sup>-3</sup></td><td>TC<sub>Z</sub>: 9.65 &times; 10<sup>-1</sup><br>TC<sub>N</sub>: 1.40 &times; 10<sup>-3</sup><br>TC<sub>Z'</sub>: 1.43 &times; 10<sup>-2</sup></td><td>1.74 &times; 10<sup>-2</sup></td><td>4.30 &times; 10<sup>-6</sup></td></tr></tbody></table>

We let f<sub>s</sub> and P<sub>s</sub> (s = 3, 2, 1, 0) represent the fractions and specific permeabilities of the +1, 0, -1 and -2 charged macrostates of TC, and f<sub>N</sub>, f<sub>Z</sub> and f<sub>Z'</sub> denote the fractions of the TC<sub>N</sub>, TC<sub>Z</sub>, and TC<sub>Z'</sub> tautomers, respectively. It naturally follows that &Sigma;<sub>s</sub> f<sub>s</sub> = 1 and f<sub>2</sub> = f<sub>N</sub> + f<sub>Z</sub> + f<sub>Z'</sub>. Using the specific permeability of TC<sub>Z</sub> as a proxy for all charged macrostates and the TC<sub>Z</sub>' tautomer, we obtain the effective permeability of TC from the pH partitioning protocol as

$$P_{\text{eff}} = \sum_{s} f_{s} P_{s} \approx (f_{0} + f_{1} + f_{Z} + f_{Z'} + f_{3}) P_{Z} + f_{N} P_{N} = (1 - f_{N}) P_{Z} + f_{N} P_{N}. \quad (4)$$

To determine f<sub>N</sub>, the ratios between the TC<sub>N</sub>, TC<sub>Z</sub> and TC<sub>Z'</sub> tautomers are needed. These pH-independent ratios cannot be determined from the macroscopic pK<sub>a</sub>s and instead require the microscopic acid dissociation constants, the negative log values of which are known as the microscopic pK<sub>a</sub>s (Fig. 2). Although their experimental measurement is challenging, estimation can be made based on model compounds that are structural analogs of selected tautomers. For instance, since tetracycline methiodide (TCMI) hosts a permanent positive charge on its quaternary ammonium at site B,<sup>66</sup> once it loses a proton on its tricarbonyl group A, the compound can be taken as an analog of A<sup>-</sup>B<sup>+</sup>C<sup>0</sup>. The second macroscopic pKa of TCMI has therefore been used by Leeson et al. as an estimate of k<sub>13</sub>, i.e., pk<sub>13</sub> ≈ 7.80.<sup>64</sup> Similarly, k<sub>21</sub> and k<sub>12</sub> have been estimated using another
model compound as well as TC's macroscopic $pK_a$s: $pk_{21} \approx 5.97$ and $pk_{12} \approx 8.8$.⁶⁴ Assuming that the acidity of the tricarbonyl group A is approximately the same in the neutral $A^0B^0C^0$ and the zwitterionic $A^0B^+C^-$, i.e., $pk_{31} \approx pk_{21}$, the ratio between the three tautomers can be readily estimated from thermodynamic cycles in the ionization scheme⁶³—one such thermodynamic cycle involving $A^-B^+C^0$, $A^-B^0C^0$ and $A^0B^0C^0$ is highlighted in bold in Fig. 2, from which $f_N/f_Z$ is found to be $k_{12}/k_{21} = 10^{5.97-8.8} = 10^{-2.83}$. Similarly, $f_{Z'}/f_Z$ is obtained from the thermodynamic cycle involving $A^-B^+C^0$, $A^-B^+C^-$ and $A^0B^+C^-$ as $f_{Z'}/f_Z = k_{13}/k_{31} = 10^{5.97-7.8} = 10^{-1.83}$. The fractions of the three tautomers TC$_N$, TC$_Z$ and TC$_{Z'}$ can then be calculated from their sum ($f_2$). These values, in particular $f_N$, are used to compute the effective permeability of TC via Eq. 4.

<table><caption>**Table 2: Specific permeabilities of TC$_N$ ($P_N$), TC$_Z$ ($P_Z$) and the effective permeability ($P_{eff}$) of TC at pH = 6 (unit: cm/s).**</caption><thead><tr><th>P<sub>N</sub></th><th>P<sub>Z</sub></th><th>P<sub>eff</sub> (pH-partitioning)</th><th>P<sub>eff</sub> (BWAP)</th></tr></thead><tbody><tr><td>2.22 &times; 10<sup>-2</sup></td><td>1.63 &times; 10<sup>-8</sup></td><td>3.17 &times; 10<sup>-5</sup></td><td>3.17 &times; 10<sup>-5</sup></td></tr></tbody></table>

At pH = 6, the above calculation yields $P_{eff} = 3.17 \times 10^{-5}$ cm/s ($\log P_{eff} = -4.50$), which is similar to the result of a PAMPA assay using a decane membrane with 10% soybean lecithin at pH = 6 ($P_{PAMPA} \approx 1.8 \times 10^{-5}$ cm/s from Fig. 3 of Yamauchi and Sugano³⁴), but differs by approximately one log unit from $\log P_{app} = -5.86$ measured using fluorescence assays conducted on POPC liposomes.²⁵,²⁶ Before we further discuss the nature of this deviation, we repeat the above calculation using the BWAP protocol: The PMF and diffusivity profiles of TC's charged macrostates and the TC$_Z$ tautomer are again approximated by those of the TC$_Z$ tautomer, yielding a weighted average potential $w_m(z)$ as:

$$e^{-\beta w_m(z)} = \sum_s f_s e^{-\beta w_s(z)} \approx (1 - f_N) e^{-\beta w_Z(z)} + f_N e^{-\beta w_N(z)}, \quad (5)$$

where $w_Z(z)$ and $w_N(z)$ represent the original PMFs of TC$_Z$ and TC$_N$, respectively (Fig. 2 A). The PMF of TC$_N$ can be shifted by $\Delta G = -1/\beta \ln(f_N/f_Z)$ to account for the free-energy cost for converting a TC$_Z$ tautomer to TC$_N$ in bulk solution at the given pH. The resulting shifted PMF of TC$_N$ enables the direct comparison of its permeation energetics with that of TC$_Z$, after the different
fractions of the two tautomers are taken into consideration. As given in Table 2, calculation results based on the BWAP and the pH-partitioning schemes are indistinguishable. This behavior arises from the significantly higher free energy barrier of $TC_Z$ than that of the shifted PMF of $TC_N$, i.e., even after the lower fraction of the neutral tautomer has been taken into account, its contribution to $P_{\text{eff}}$ remains dominant. As pointed out by Harris et al.,¹⁸ that when one ionization state dominates $P_{\text{eff}}$, both approaches yield similar numerical results.

The effective permeability of TC at pH = 6 depends minimally on the specific permeabilities of its charged macrostates of TC and that of the $TC_Z$ tautomer: setting these specific permeabilities to zero yields negligible reduction in $P_{\text{eff}}$; similarly, increasing them by three orders of magnitude raises the computed $P_{\text{eff}}$ by merely ~1.7% (Table S2). Only when their specific permeabilities become five orders of magnitude greater, i.e., approaching that of $P_N$, the computed $P_{\text{eff}}$ increases by approximately 2.7 fold (Table S2). However, as its zwitterionic components have rendered the specific permeability of $TC_Z$ much lower than its neutral tautomer, neither $TC_Z$ nor any charged macrostate of TC is expected to have specific permeability approaching that of $TC_N$. The overall insensitivity of the computed $P_{\text{eff}}$ to specific permeabilities of the charged macrostates is explained by the dominant fraction of the zero-net-charge macrostate at pH = 6, as well as the large specific permeability of $TC_N$. Since $P_N$ is six orders of magnitude greater than $P_Z$, the neutral $TC_N$ dominates $P_{\text{eff}}$, even though its fraction is only ~1/1000 that of $TC_Z$. Indeed, as shown in Table S3, $TC_N$ makes the dominant contribution to $P_{\text{eff}}$ in acidic to basic environments with pH ranging from 1 to 10. However, we should add that this conclusion relies on the accuracy of the estimated $f_N/f_Z$ ratio, which is based on microscopic acid dissociation constants inferred from analog compounds. As discussed further below, the choice of these model compounds and the uncertainty in their macroscopic $pK_a$ values directly affect the reliability of $f_N/f_Z$, and, hence, of the computed membrane permeability.
## Bilayer patch size significantly affects TC<sub>N</sub> permeation energetics

While the TC<sub>N</sub> tautomer faces a rather moderate energy barrier against its permeation, trajectory analysis reveals significant bilayer deformation as the molecule traverses the membrane. As shown in Fig. S3, the monolayer occupied by TC<sub>N</sub> “caves in” near the bilayer center, with deformation spanning over half the (x, y) plane. This observation motivates an investigation into the impact of bilayer patch size on TC permeability estimates, i.e., whether a minimum patch size is required to accommodate TC-induced membrane deformation and provide unbiased permeation energetics. Since TC<sub>N</sub> dominates the effective permeability of the drug molecule over a wide range of pH, we focus our analysis on this tautomer. Towards this end, we performed additional PMF calculations for TC<sub>N</sub> in bilayers containing 32, 64, and 256 POPC molecules. Together with the previously discussed 128-POPC system, these simulations enable a systematic assessment of patch-size effects. For comparison, we also simulated three other permeants—methanol, urea, and DDA—that differ in size and polarity (Table S4).<sup>67</sup> The convergence of these PMF calculations is illustrated in Fig. S4 and Fig. S5.

As shown in Fig. 5 and Table S5, patch size exerts minimal influence on methanol, urea, and DDA, with permeation barriers differing by at most ~9% across patch sizes. Specifically, methanol has a barrier of 3.4 kcal/mol in the 32-POPC system and 3.6 kcal/mol in the remaining three lipid patches—consistent with previous work reporting 3.5 ± 0.2 kcal/mol in a 100-POPC bilayer at 308 K.<sup>58</sup> Urea exhibits no significant energy well, with a barrier ranging from 10.0 to 10.6 kcal/mol across the four systems, in line with the 9.4 kcal/mol barrier reported for a DMPC bilayer at 298.15 K.<sup>6</sup> As the patch size increases, the barrier region ($w(z) > 0$) becomes slightly narrower—a trend also reported by Lee et al.<sup>6</sup> Finally, DDA shows barriers between 7.4 and 8.1 kcal/mol, comparable to the 7.5 kcal/mol value observed for a 128-POPC bilayer at 308 K.<sup>68</sup>

By contrast, the permeation energetics of TC<sub>N</sub> exhibits a pronounced patch-size dependence. The energy barrier decreases from 8.1 kcal/mol in the 32-POPC patch to 7.0 kcal/mol in the 64-POPC patch, 4.4 kcal/mol in the 128-POPC system, and just 2.2 kcal/mol in the 256-POPC system. To understand this dependence, we compare the hydrogen-bonding interactions, hydration levels,
<figure><img src="image_13.png" /><figcaption>Figure 5: Membrane permeation PMFs and analysis of hydrogen-bond interactions and the number of solvating water molecules for TC<sub>N</sub>, methanol, urea, and DDA. (A) Symmetrized PMFs for each permeant. (B) Number of hydrogen bonds between each permeant and POPC. (C) Number of hydrogen bonds between each permeant and water. (D) Number of water molecules within 2.4 Å of each permeant. The error bars represent the standard deviations.</figcaption></figure>
and bilayer deformation across the four systems. TCN forms on average 1.5 hydrogen bonds with POPC at the lipid-water interface (z = ±15 Å), exceeding values for methanol (0.2), urea (1.0), and DDA (0.6). These stronger interactions help explain why TCN, despite its size and Ro5 violation, has a lower barrier than urea or DDA. To minimize the energetic cost of membrane crossing, TCN recruits hydrogen bonding partners from lipids or water, with the success of this recruitment varying with patch size. For example, at z = 0, TCN maintains about 1 hydrogen bond in 128- and 256-POPC systems, but only about 0.5 in the 64-POPC patch, and virtually none in the 32-POPC system. Thermal undulations, which support such interactions, are constrained by simulation box size,⁶⁹ making smaller bilayers effectively stiffer. Formation of hydrogen bonds at z = 0 demands substantial membrane deformation (Fig. 4) and increased solvation deep in the bilayer (Fig. 5). Smaller patches, like 32-POPC, resist such deformation, producing artificially high barriers. As shown in Fig. S6 and Fig. S7, membrane thickness reduction is minimal for the 32-POPC bilayer, but evident in the larger systems. The inadequacy of the smallest 32-POPC patch in accommodating TCN is further illustrated in the flip-flop behavior of the permeant. While an average orientation of ~90° is recorded at z = 0 in all four bilayer patches (Fig. S8), the distributions of the $\vec{v}_2$ tilt angle differ significantly: In both the 128-POPC and 256-POPC patches, and, to a lesser degree, the 64-POPC patch, TCN mostly maintains a “horizontal orientation” where the angle between $\vec{v}_2$ and the z-axis is around 30° or 150° even in the bilayer center region (−1 < z < 1 Å); in sharp contrast, the same tilt angle is ~90° in the 32-POPC bilayer (Fig. S9), indicating that the smallest patch cannot support TCN in its energetically favorable orientation. Along with the hydrogen-bonding patterns and membrane thickness profiles, this result explains the elevated energy barrier observed in the smaller bilayer systems for TCN.

While the smaller bilayers present artificially high barriers against TCN permeation, further analysis suggests that the larger ones can be plagued by hysteresis. This issue is particularly evident once the patch size increases to 256 POPC. As illustrated in Fig. 6, a TCN molecule traversing the membrane from the lower to the upper leaflet may not undergo its flip-flop transition even at $z \approx 11.1$ Å in the 256-POPC bilayer. The phenolic diketone side of TCN continues to be solvated by
<figure><img src="image_14.png" /><figcaption>Figure 6: Snapshots of TC~N~ in the 256-POPC system. (A) TC~N~ entering the bilayer from the upper leaflet. (B) TC~N~ entering the bilayer from the lower leaflet. The location of TC~N~ and the angle between its $\vec{v}_2$ vector and the z-axis are 13.7 Å and 41° in (A), 11.1 Å and 138° (t = 28.2 ns), 15.0 Å and 99° (t = 28.3 ns), as well as 13.8 Å and 28° (t = 38.2 ns) in (B), respectively. Water is rendered as semi-transparent surfaces, with those near TC~N~ highlighted in van der Waals representation. POPC molecules are shown in semi-transparent licorice style, with their phosphorus atoms shown as tan spheres.</figcaption></figure>

water molecules that enter from the lower leaflet side. Its flip-flop transition only occurs at approximately $z = 15$ Å, after which the $\vec{v}_2$ tilt angle eventually switches to ~30°, as typically observed at this location (Fig. 6 B). Such a lag in the change of TC~N~ orientation relative to its z-coordinate is a clear sign of hysteresis. As conformations where the membrane is minimally deformed and those with severe membrane deformation (Fig. 6) all contribute to the force samples collected at $z \approx 11-15$ Å, without proper weighting this z-range may be rendered less favorable than it actually is. Indeed, while the other three bilayers consistently reveal an energy well of approximately -4.9 kcal/mol at $z = \pm 14$ Å, the 256-POPC system shows a shallower well of -4.3 kcal/mol around this location. These values suggest that the impact of sampling hysteresis is the most severe in the largest bilayer patch. As noted earlier, the 256-POPC bilayer is the softest among the four, which explains its greatest hysteresis effect observed above. In theory, such hysteresis can be resolved by increasing simulation time to improve the estimate of the relative contribution from the two types of conformations shown in Fig. 6 A and Fig. 6 B. In practice, the large system size of the 256-POPC bilayer limits its sampling efficiency—apart from an increased computational cost, the long-wavelength membrane deformations coupled with the permeant-lipid interactions suggests that multi-dimensional PMF calculations, where the permeant orientation and/or even the membrane deformation serve as additional collective variable(s), may be needed. While such
calculations are beyond the scope of the current study, it appears that for one-dimensional PMF calculations, the 128-POPC patch has highly similar membrane thickness profile and hydrogen bonding pattern to the 256-POPC (Fig. S6 and Fig. 5). Therefore, the former appears to adequately accommodate TC<sub>N</sub>-induced membrane deformation without significant hysteresis observed in the latter bilayer. Taken together, our results suggest that while a membrane that is too small should be avoided, e.g., the 32-POPC bilayer, lipid patches are not necessarily “the larger, the better”. A bilayer that is large enough to accommodate deformation induced by the permeant without introducing severe sampling hysteresis should be employed in its one-dimensional PMF calculations.

Discussion

TC is one of the most widely used broad-spectrum antibiotics.²⁹ Administered orally, it acts by inhibiting protein synthesis through targeting the bacterial ribosome. The amount of absorbed TC is directly proportional to its administered dose over a tenfold range,³² indicating that TC absorption is a passive process.³⁵ In both PAMPA and liposomal fluorescence assays, TC demonstrates moderate permeability on the order of 10⁻⁶ to 10⁻⁵ cm/s.²⁵,²⁶,³⁴ These values are consistent with its rapid adsorption and good *in vivo* bioavailability,³²,³³ but appear paradoxical, given the six hydrogen-bond donors of the drug and its dominant zwitterionic population at neutral pH. Through permeation energetics and diffusivity calculations, we show that unlike the previously hypothesized cancellation-of-charge or ion-pair mechanism³⁴,³⁵, the unexpectedly decent permeability of TC is attributed to its neutral tautomer TC<sub>N</sub>. While it has a fraction only ~1/1000 that of the zwitterionic TC<sub>Z</sub>, its specific permeability is approximately six orders of magnitude greater. Consequently, the effective permeability of TC is dominated by TC<sub>N</sub> over a wide range of pH values (Table S3).

To calculate the effective permeability of TC, one must combine the contributions of its various protonation and tautomeric states using protocols such as pH-partitioning or BWAP.¹⁸ Although they differ in their assumptions regarding how fast protonation is relative to membrane permeation, the two protocols yield indistinguishable numerical results for TC. This behavior arises from the
significant difference between the PMF of $TC_Z$ and the shifted PMF of $TC_N$. As shown in Fig.2, even after their population fractions are taken into account, the permeation barrier faced by $TC_N$, as revealed by its shifted PMF, is still much lower than $TC_Z$. Therefore, whether the PMFs are first integrated to produce the specific permeabilities and then combined to yield the effective permeability as in the pH-partitioning protocol, or, combined first to yield an effective PMF and then used to compute $P_{eff}$ as in the BWAP treatment, the numerical results are effectively indistinguishable. In other words, the population of the $TC_Z$ tautomers near the bilayer center is sufficiently low such that their membrane permeation behavior—whether they cross as zwitterionic $TC_Z$ or switch to neutral $TC_N$—becomes insignificant in determining the value of $P_{eff}$. Despite this numerical result, it is still of interest to ask, which one of these two processes, protonation or membrane permeation, is faster for TC? While an extensive analysis based on the transition path theory as carried out by Harris et al.¹⁸ is beyond the scope of this work, we note that the transition between $TC_Z$ and $TC_N$ does not involve any net release or uptake of protons. The switch between the two tautomers can proceed by transferring a proton from one ionizable group to another, i.e., the molecule serves simultaneously as its own proton donor and acceptor. Therefore, the low concentration of water within the membrane, which can slow down the protonation or deprotonation that involves net proton exchange, does not hinder this tautomeric transition. Given that TC only has a moderate permeability, it is plausible that the transition between its zwitterionic and neutral tautomers occurs faster than its membrane permeation.

As pointed out by Sezer and Oruç,²² results of a liposomal fluorescence assay may become unreliable if, instead of membrane permeation, the protonation/deprotonation of the permeant or the fluorescent dye becomes the rate-limiting step. At pH = 6, TC predominantly carries a zero net charge, and its macrostate carrying a −1 charge corresponds to the second highest fraction. In the fluorescence assay detailed previously,⁸,²⁶ the molecule acts as a weak acid upon entering the liposome, i.e., TC crosses the membrane carrying zero net charge and then releases a proton. The pH-sensitive dye HPTS enclosed in the liposome then captures a proton, which in turn reduces its fluorescence signal. The rate of the fluorescence reduction ($k$) is used to estimate the apparent
permeability of TC according to $P_{\text{app}} = kd/6$, where $d \approx 167.7$ nm is the mean hydrodynamic diameter of the POPC liposomes.²⁶ Following Harris et al.,¹⁸ we estimate the on-rate of protonation by $k_{\text{on}} = k_0 \times 10^{-\text{pH}}$, where $k_0 = 10^{10} \text{ M}^{-1}\text{s}^{-1}$ is the intrinsic protonation rate for small molecules. With the liposomal fluorescence measurement performed at pH = 6, we have $k_{\text{on}} = 10^{10-6} = 10^4 \text{ s}^{-1}$, which is taken to be the protonation rate for both HPTS and TC in the pH = 6 liposomal environment. The deprotonation rate $k_{\text{off}}$ is estimated from $k_{\text{on}}$ and the corresponding acid dissociation constant.¹⁸ Specifically, we use $k_{\text{off}} = k_0 \times 10^{-\text{p}K_{\text{a2}}}$ to estimate the rate of deprotonation of the zero-net-charge TC, obtaining $k_{\text{off}} = 10^{10-7.75} = 10^{2.25} \text{ s}^{-1}$. In comparison, the measured apparent permeability of TC ($P_{\text{app}} = 10^{-5.86} \text{ cm/s}$) corresponds to an intrinsic rate $k = 6P_{\text{app}}/d \approx 0.5 \text{ s}^{-1}$. Since both $k_{\text{on}}$ and $k_{\text{off}}$ exceed $k$ by orders of magnitude, the rate of fluorescence change cannot be reflecting that of TC deprotonation or HPTS protonation. Thus, unlike propranolol,²² the liposomal fluorescence assay result for TC should reflect the permeation dynamics of the drug molecule.

The many hydrogen bond donors and acceptors of TC present some unique challenges to its permeability calculation. For instance, several of its six hydrogen bond donors are part of the two ionizable groups of the molecule. With altogether three ionizable groups, the complete ionization scheme of TC involves eight microstates, which belong to four macrostates distinguished by the net charge carried by the molecule. Calculation of TC's effective permeability requires not only the fraction of each macrostate, which can be readily computed using its macroscopic $pK_{\text{a}}$, but also the ratios between the tautomers that correspond to microstates in a given macrostate. These ratios depend on microscopic $pK_{\text{a}}$, which are considerably more difficult to estimate than macroscopic ones. While assignment based on analogous model compounds is a well-established strategy, the uncertainty associated to the $pK_{\text{a}}$ measurement of these model compounds, and/or the choice of different model compounds may both affect the computed $P_{\text{eff}}$. As an example, had $pK_{21}$ been estimated from a different, slightly more dissimilar compound,⁷⁰,⁷¹ an up to 0.5 log unit reduction in $P_{\text{eff}}$ would be produced (Table S6).

Compared to the uncertainty in microscopic $pK_{\text{a}}$ estimation, incorrect use of macroscopic $pK_{\text{a}}$
pK<sub>a1</sub> ≈ 3.3 to the tricarbonyl group A, conflicting assignments have been made for the remaining two pK<sub>a1</sub>s: Stephens et al.<sup>31</sup> assigned pK<sub>a2</sub> to the dimethylammonium group B and pK<sub>a3</sub> to the phenolic diketone group C, while the opposite assignment was given by Leeson et al.<sup>64</sup> More recently, Qiang and Adams have shown that the latter two ionizable groups undergo deprotonation simultaneously during titration.<sup>65</sup> As a result, pK<sub>a2</sub> and pK<sub>a3</sub> both have significant contributions from the dimethylammonium and the phenolic diketone groups, and neither of the two macroscopic pK<sub>a1</sub>s can be unambiguously assigned to either group. Naively assigning the macroscopic K<sub>a2</sub> and K<sub>a3</sub> to the microscopic k<sub>2</sub> and k<sub>3</sub> (or in the reverse order) changes the computed P<sub>eff</sub> by up to three orders of magnitude (Table S7), reflecting the substantial error introduced by incorrect macroscopic pK<sub>a</sub> assignments. While the ionization scheme of TC appears more complex than a simple acid, other drug molecules, especially the bRo5 ones, may well present similarly (if not more) complex ionization schemes. Therefore, care should be taken to avoid the incorrect use of macroscopic pK<sub>a1</sub>s in permeability calculation of these molecules.

Another challenge presented by TC is its strong hydrogen bonding capacity with lipid molecules. As it elicits surrounding lipids to help coordinating its many hydrogen bond donors and acceptors, membranes that cannot readily deform and cater to TC's hydrogen bonding requests due to finite-size effects may artificially increase energy barriers against its permeation. As noted earlier, this result does not mean that larger membranes should always be preferred in the permeability calculation of TC and other drugs with strong hydrogen bonding capacity. As demonstrated by the 256-POPC bilayer, significant hysteresis may occur in calculations using this large membrane. Altogether, our results indicate that both the smallest 32-POPC and the largest 256-POPC bilayers should be avoided here. Between the remaining two patches, the 128-POPC bilayer appears to better accommodate deformation induced by TC without introducing severe sampling hysteresis. Nonetheless, further studies with a broader range of patch sizes, sampling protocols, and/or multidimensional PMF calculations—or, alternatively, a path search in multidimensional CV subspace
combined with a one-dimensional PMF calculation along the discovered path⁷²—will be necessary
to determine the onset of hysteresis and explore efficient mitigation strategies.

Conclusion

Molecular simulations have long been employed to estimate the membrane permeability and probe
the passive transport mechanism of drug or drug-like small molecules. The simulated drug molecules
typically obey Lipinski's "Rule-of-Fve", a set of guidelines for screening compounds with desirable bioavailability. Whether these simulations can tackle more complex permeants, especially
those beyond-Rule-of-Five drugs that violate more than one Ro5, remains largely unknown. With
six hydrogen bond donors, the antibiotic TC violates one rule under the Ro5. While this still
places it within the class of orally active drugs having "no more than one violation of the Ro5",
TC represents an ideal borderline case to explore the capacity of molecular simulations in tackling
more complex and challenging permeants. In this work, we calculated the PMFs and diffusivity
profiles governing the membrane permeation of two zero-net-charge tautomers of TC—the neutral
TC<sub>N</sub> and the zwitterionic TC<sub>Z</sub>—across 128-POPC lipid bilayers, and determined their intrinsic,
pH-independent specific permeabilities using the ISD model. While the diffusivity profiles of TC<sub>N</sub>
and TC<sub>Z</sub> are highly similar, their PMFs differ substantially, with the free-energy barrier of TC<sub>N</sub>
being 9 kcal/mol lower than that of TC<sub>Z</sub>. This large difference leads to a six-order-of-magnitude
higher specific permeability of TC<sub>N</sub> than TC<sub>Z</sub>. As a result, despite a fraction only ~1/1000 that of
the zwitterionic tautomer, the neutral TC<sub>N</sub> dominates the effective permeability of TC over a wide
pH range.

The computed effective permeability of TC at pH = 6 is 3.17 × 10⁻⁵ cm/s, which is close to the
result of a PAMPA assay (1.8 × 10⁻⁵ cm/s) performed on a decane membrane with 10% soybean
lecithin,³⁴ but differs by approximately one log unit from that measured using a liposomal fluorescence assay (log P<sub>app</sub> = −5.86).²⁵,²⁶ Overall, the approximately one log unit deviation between the
computationally estimated P<sub>eff</sub> and the experimentally measured P<sub>app</sub> is comparable to the perfor-
mance of permeability calculations for other drug and drug-like molecules.<sup>9,11,17,20</sup> However, this result depends closely on the pH-independent ratio between $TC_N$ and $TC_Z$ ($f_N/f_Z$), which is estimated from the corresponding microscopic acid dissociation constants. Our calculations indicate that when incorrectly assigned macroscopic $pK_{as}$ are used to estimate $f_N/f_Z$, $P_{eff}$ may differ by orders of magnitude from that obtained using judiciously estimated microscopic $pK_{as}$. As ionization schemes can become particular complex in drugs violating the Ro5 on the number of hydrogen bond donors and/or acceptors, extra care should be taken to avoid the incorrect use of macroscopic $pK_{as}$ in their permeability calculations.

The observation of significant membrane deformation—up to half of the simulation box in the $(x, y)$ plane—prompted a systematic investigation of patch-size effects on the permeation energetics of $TC_N$, given its dominant contribution to $P_{eff}$. To provide a comparative perspective, we also examined three other molecules—methanol, urea, and DDA—which differ in size and physicochemical properties. PMFs were generated for bilayers composed of 32, 64, 128, and 256 POPC lipids and the results reveal a remarkable sensitivity of the energy barrier for $TC_N$ to patch size, in stark contrast to the generally consistent PMFs observed for the other three permeants. Specifically, the free-energy barrier for $TC_N$ decreases with increasing patch size, a trend attributable to the enhanced hydrogen bonding between $TC_N$ and lipids facilitated by the greater membrane deformation achievable in larger bilayers. Such deformation in the largest 256-POPC bilayer produces significant hysteresis that compromises the quality of the free-energy calculation. As a result, both the smallest 32-POPC and, counterintuitively, the largest 256-POPC bilayers should be avoided for $TC_N$ permeability calculation. In general, for bRo5 drugs with strong hydrogen bonding capacity, bilayers that are large enough to accommodate deformation induced by the permeants without introducing severe sampling hysteresis should be employed.

The above findings prompt broader reflection on the applicability and limitations of molecular simulation approaches to the study of complex permeants. While the agreement between simulation and experiment for the effective membrane permeability to TC is encouraging, the heavy reliance on a rare neutral species—one accounting for just a small fraction of the equilibrium
population—raises questions about the generality and robustness of the mechanism. Does the observed behavior reflect a universal principle applicable to other complex, zwitterionic drugs, or is it an idiosyncrasy of the tautomeric landscape of TC? The pronounced sensitivity of the PMF to membrane patch size for TCN also underscores the challenge of simulating the mechanically flexible bilayers, especially for large, polarizable permeants capable of remodeling their environment. These effects are rarely examined systematically, but may critically impact predictions of membrane permeabilities to drug molecules from first principles. As computational methods continue to mature, it will be essential to recognize that achieving numerical agreement with experimental observables does not guarantee mechanistic accuracy—particularly when conclusions hinge on rare or extreme configurations.

## Acknowledgement

Y.W. is supported by projects 14323816 and AoE/P-705/16 from the Hong Kong Research Grants Council as well as direct grants from the Chinese University of Hong Kong. C.C. is indebted to the CNRS for an International Emerging Action with Hong Kong.

## References

(1) Hannesschlaeger, C.; Horner, A.; Pohl, P. Intrinsic membrane permeability to small molecules. *Chemical reviews* **2019**, *119*, 5922–5953.

(2) Lipinski, C. A.; Lombardo, F.; Dominy, B. W.; Feeney, P. J. Experimental and Computational Approaches to Estimate Solubility and Permeability in Drug Discovery and Development Settings. *Advanced Drug Delivery Reviews* **2001**, *46*, 3–26.

(3) Kansy, M.; Senner, F.; Gubernator, K. Physicochemical High Throughput Screening: Parallel Artificial Membrane Permeation Assay in the Description of Passive Absorption Processes. *Journal of Medicinal Chemistry* **1998**, *41*, 1007–1010.
(4) Artursson, P.; Palm, K.; Luthman, K. Caco-2 Monolayers in Experimental and Theoretical Predictions of Drug Transport. *Advanced Drug Delivery Reviews* **2001**, *46*, 27–43.

(5) Hubatsch, I.; Ragnarsson, E. G. E.; Artursson, P. Determination of Drug Permeability and Prediction of Drug Absorption in Caco-2 Monolayers. *Nature Protocols* **2007**, *2*, 2111–2119.

(6) Christopher T. Lee; Lee, C. T.; Jeffrey Comer; Comer, J.; Conner Herndon; Herndon, C.; Nelson Leung; Nelson W.Y. Leung; Leung, N.; Ànna Pavlova et al. Simulation-Based Approaches for Determining Membrane Permeability of Small Compounds. *Journal of Chemical Information and Modeling* **2016**, *56*, 721–733.

(7) Menichetti, R.; Kanekal, K. H.; Bereau, T. Drug-Membrane permeability across chemical space. *ACS Central Sci.* **2019**, *5*, 290–298.

(8) Krämer, A.; Ghysels, A.; Wang, E.; Venable, R. M.; Klauda, J. B.; Brooks, B. R.; Pastor, R. W. Membrane permeability of small molecules from unbiased molecular dynamics simulations. *J. Chem. Phys.* **2020**, *153*, 124107.

(9) Richard M. Venable; Venable, R. M.; Andreas Krämer; Krämer, A.; Richard W. Pastor; Pastor, R. W. Molecular Dynamics Simulations of Membrane Permeability. *Chemical Reviews* **2019**, *119*, 5954–5997, MAG ID: 2911593691.

(10) 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.

(11) Harris, J.; Roux, B. Membrane Permeability of Sucrose Calculated from Equilibrium Time-Correlation Functions Using Molecular Dynamics Simulations with Enhanced Sampling.

(12) Diamond, J. M.; Katz, Y. Interpretation of Nonelectrolyte Partition Coefficients between Dimyristoyl Lecithin and Water. *The Journal of Membrane Biology* **1974**, *17*, 121–154.
(13) Marrink, S.-J.; Berendsen, H. J. C. Simulation of Water Transport through a Lipid Membrane.
*The Journal of Physical Chemistry* **1994**, *98*, 4155–4168.

(14) Sun, R.; Dama, J. F.; Tan, J. S.; Rose, J. P.; Voth, G. A. Transition-Tempered Metadynamics Is a Promising Tool for Studying the Permeation of Drug-like Molecules through Membranes.
*Journal of Chemical Theory and Computation* **2016**, *12*, 5157–5169.

(15) Tse, C. H.; Comer, J.; Chu, S. K. S.; Wang, Y.; Chipot, C. Affordable Membrane Permeability Calculations: Permeation of Short-Chain Alcohols through Pure-Lipid Bilayers and a Mammalian Cell Membrane.
*J. Chem. Theory Comput.* **2019**, *15*, 2913–2924.

(16) Shoji, A.; Kang, C.; Fujioka, K.; Rose, J. P.; Sun, R. Assessing the Intestinal Permeability of Small Molecule Drugs via Diffusion Motion on a Multidimensional Free Energy Surface.
*Journal of Chemical Theory and Computation* **2022**, *18*, 503–515.

(17) Chipot, C. Predictions from First-Principles of Membrane Permeability to Small Molecules: How Useful Are They in Practice?
*Journal of Chemical Information and Modeling* **2023**, *63*, 4533–4544.

(18) Harris, J.; Chipot, C.; Roux, B. How Is Membrane Permeation of Small Ionizable Molecules Affected by Protonation Kinetics?
*The Journal of Physical Chemistry B* **2024**, *128*, 795–811.

(19) Kang, C.; Shoji, A.; Chipot, C.; Sun, R. Impact of the Unstirred Water Layer on the Permeation of Small-Molecule Drugs.
*Journal of Chemical Information and Modeling* **2024**, *64*, 933–943.

(20) Tse, C. H.; Comer, J.; Wang, Y.; Chipot, C. Link between Membrane Composition and Permeability to Drugs.
*J. Chem. Theory Comput.* **2018**, *14*, 2895–2909.

(21) Yue, Z.; Li, C.; Voth, G. A.; Swanson, J. M. Dynamic Protonation Dramatically Affects the Membrane Permeability of Drug-like Molecules.
*J. Am. Chem. Soc.* **2019**, *141*, 13421–13433.
(22) Sezer, D.; Oruç, T. Protonation Kinetics Compromise Liposomal Fluorescence Assay of Membrane Permeation. *The Journal of Physical Chemistry B* **2017**, *121*, 5218–5227.

(23) Eyer, K.; Paech, F.; Schuler, F.; Kuhn, P.; Kissner, R.; Belli, S.; Dittrich, P. S.; Krämer, 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.

(24) Barry, P. H.; Diamond, J. M. Effects of unstirred layers on membrane phenomena. *Physiol. Rev.* **1984**, *64*, 763–872.

(25) Krämer, S. D.; Aschmann, H. E.; Hatibovic, M.; Hermann, K. F.; Neuhaus, C. S.; Brunner, C.; Belli, S. When barriers ignore the “rule-of-five”. *Adv. Drug Deliv. Rev.* **2016**, *101*, 62–74.

(26) Hermann, K. F.; Neuhaus, C. S.; Micallef, V.; Wagner, B.; Hatibovic, M.; Aschmann, H. E.; Paech, F.; Alvarez-Sanchez, R.; Krämer, S. D.; Belli, S. Kinetics of Lipid Bilayer Permeation of a Series of Ionisable Drugs and Their Correlation with Human Transporter-Independent Intestinal Permeability. *European Journal of Pharmaceutical Sciences: Official Journal of the European Federation for Pharmaceutical Sciences* **2017**, *104*, 150–161.

(27) Lipinski, C. A.; Lombardo, F.; Dominy, B. W.; Feeney, P. J. Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. *Adv. Drug Deliv. Rev.* **1997**, *23*, 3–25.

(28) Hartung, I. V.; Huck, B. R.; Crespo, A. Rules were made to be broken. *Nature Reviews Chemistry* **2023**, *7*, 3–4.

(29) Chopra, I.; Roberts, M. Tetracycline antibiotics: mode of action, applications, molecular biology, and epidemiology of bacterial resistance. *Microbiology and molecular biology reviews* **2001**, *65*, 232–260.

(30) Zhang, C.; Ouyang, Q.; Zhou, X.; Huang, Y.; Zeng, Y.; Deng, L.; Lin, D.; Zheng, W. In
vitro activity of tetracycline analogs against multidrug-resistant and extensive drug resistance clinical isolates of Mycobacterium tuberculosis. *Tuberculosis* **2023**, *140*, 102336.

(31) Stephens, C. R.; Murai, K.; Brunings, K. J.; Woodward, R. B. Acidity Constants of the Tetracycline Antibiotics. *J. Am. Chem. Soc.* **1956**, *78*, 4155–4158.

(32) Pindell, M.; Cull, K.; Doran, K.; Dickison, H. Absorption and excretion studies on tetracycline. *The Journal of Pharmacology and Experimental Therapeutics* **1959**, *125*, 287–294.

(33) Tam, K. Y.; Avdeef, A.; Tsinman, O.; Sun, N. The Permeation of Amphoteric Drugs through Artificial Membranes–an in Combo Absorption Model Based on Paracellular and Transmembrane Permeability. *Journal of Medicinal Chemistry* **2010**, *53*, 392–401.

(34) Yamauchi, S.; Sugano, K. Permeation characteristics of tetracyclines in parallel artificial membrane permeation assay. *ADMET and DMPK* **2019**, *7*, 151–160.

(35) Colaizzi, J. L.; Klink, P. R. pH-Partition Behavior of Tetracyclines. *Journal of Pharmaceutical Sciences* **1969**, *58*, 1184–1189.

(36) Venable, R. M.; Brown, F. L.; Pastor, R. W. Mechanical properties of lipid bilayers from molecular dynamics simulation. *Chemistry and physics of lipids* **2015**, *192*, 60–74.

(37) Venable, R. M.; Ingólfsson, H. I.; Lerner, M. G.; Perrin, B. S.; Camley, B. A.; Marlink, S. J.; Brown, F. L. H.; Pastor, R. W. Lipid and Peptide Diffusion in Bilayers: The Saffman–Delbrück Model and Periodic Boundary Conditions. *The Journal of Physical Chemistry B* **2017**, *121*, 3443–3457.

(38) Castro-Román, F.; Benz, R. W.; White, S. H.; Tobias, D. J. Investigation of finite system-size effects in molecular dynamics simulations of lipid bilayers. *The Journal of Physical Chemistry B* **2006**, *110*, 24157–24164.

(39) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: A Web-Based Graphical User Interface for CHARMM. *Journal of Computational Chemistry* **2008**, *29*, 1859–1865.
(40) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kalé, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. *Journal of Computational Chemistry* **2005**, *26*, 1781–1802.

(41) Phillips, J. C.; Hardy, D. J.; Maia, J. D. C.; Stone, J. E.; Ribeiro, J. V.; Bernardi, R. C.; Buch, R.; Fiorin, G.; Hénin, J.; Jiang, W. et al. Scalable Molecular Dynamics on CPU and GPU Architectures with NAMD. *The Journal of Chemical Physics* **2020**, *153*, 044130.

(42) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W. Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. *J. Phys. Chem. B* **2010**, *114*, 7830–7843.

(43) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. *The Journal of Chemical Physics* **1983**, *79*, 926–935.

(44) Yu, W.; He, X.; Vanommeslaeghe, K.; MacKerell Jr, A. D. Extension of the CHARMM general force field to sulfonyl-containing compounds and its utility in biomolecular simulations. *Journal of computational chemistry* **2012**, *33*, 2451–2468.

(45) Vanommeslaeghe, K.; Raman, E. P.; MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) II: Assignment of Bonded Parameters and Partial Atomic Charges. *J. Chem. Inf. Model.* **2012**, *52*, 3155–3168.

(46) Aleksandrov, A.; Simonson, T. The Tetracycline: Mg2+ Complex: A Molecular Mechanics Force Field. *Journal of Computational Chemistry* **2006**, *27*, 1517–1533.

(47) Aleksandrov, A.; Simonson, T. Molecular Mechanics Models for Tetracycline Analogs. *J. Comput. Chem.* **2009**, *30*, 243–255.
(48) Jiang, W.; Hardy, D. J.; Phillips, J. C.; MacKerell, A. D.; Schulten, K.; Roux, B. High-Performance Scalable Molecular Dynamics Simulations of a Polarizable Force Field Based on Classical Drude Oscillators in NAMD. *J. Phys. Chem. Lett.* **2011**, *2*, 87–92.

(49) Feller, S. E.; Zhang, Y.; Pastor, R. W.; Brooks, B. R. Constant pressure molecular dynamics simulation: The Langevin piston method. *The Journal of Chemical Physics* **1995**, *103*, 4613–4621.

(50) Darden, T.; York, D.; Pedersen, L.; others Particle mesh Ewald: An N log (N) method for Ewald sums in large systems. *Journal of chemical physics* **1993**, *98*, 10089–10089.

(51) Tuckerman, M.; Berne, B. J.; Martyna, G. J. Reversible Multiple Time Scale Molecular Dynamics. *The Journal of Chemical Physics* **1992**, *97*, 1990–2001.

(52) Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. *Journal of Computational Chemistry* **1992**, *13*, 952–962.

(53) Andersen, H. C. Rattle: A “Velocity” Version of the Shake Algorithm for Molecular Dynamics Calculations. *Journal of Computational Physics* **1983**, *52*, 24–34.

(54) Fu, H.; Shao, X.; Cai, W.; Chipot, C. Taming Rugged Free Energy Landscapes Using an Average Force. *Accounts of Chemical Research* **2019**, *52*, 3254–3264.

(55) Fu, H.; Chen, H.; Wang, X.; Chai, H.; Shao, X.; Cai, W.; Chipot, C. Finding an Optimal Pathway on a Multidimensional Free-Energy Landscape. *Journal of Chemical Information and Modeling* **2020**, *60*, 5366–5374.

(56) Miao, M.; Fu, H.; Zhang, H.; Shao, X.; Chipot, C.; Cai, W. Avoiding non-equilibrium effects in adaptive biasing force calculations. *Mol. Simul.* **2021**, *47*, 390–394.

(57) Comer, J.; Chipot, C.; González-Nilo, F. D. Calculating Position-Dependent Diffusivity in Biased Molecular Dynamics Simulations. *Journal of Chemical Theory and Computation* **2013**, *9*, 876–882.
(58) Chipot, C.; Comer, J. Subdiffusion in Membrane Permeation of Small Molecules. *Scientific Reports* **2016**, *6*, 35913.

(59) Comer, J.; Schulten, K.; Chipot, C. Permeability of a Fluid Lipid Bilayer to Short-Chain Alcohols from First Principles. *J. Chem. Theory Comput.* **2017**, *13*, 2523–2532.

(60) Humphrey, W.; Dalke, A.; Schulten, K. VMD: visual molecular dynamics. *Journal of Molecular Graphics* **1996**, *14*, 33–38, 27–28.

(61) Metzler, R.; Klafter, J. The random walk's guide to anomalous diffusion: a fractional dynamics approach. *Physics reports* **2000**, *339*, 1–77.

(62) Qi, Y.; Fung, L. Y.; Chipot, C.; Wang, Y. Probing the orientation and membrane permeation of rhodamine voltage reporters through molecular simulations and free energy calculations. *Journal of Materials Chemistry B* **2025**, *13*, 2015–2028.

(63) Gunner, M. R.; Murakami, T.; Rustenburg, A. S.; Işık, M.; Chodera, J. D. Standard State Free Energies, Not pKas, Are Ideal for Describing Small Molecule Protonation and Tautomeric States. *Journal of Computer-Aided Molecular Design* **2020**, *34*, 561–573.

(64) Leeson, L. J.; Krueger, J. E.; Nash, R. A. Concerning the structural assignment of the second and third acidity constants of the tetracycline antibiotics. *Tetrahedron Lett.* **1963**, *4*, 1155–1160.

(65) Qiang, Z.; Adams, C. Potentiometric determination of acid dissociation constants (pKa) for human and veterinary antibiotics. *Water Res.* **2004**, *38*, 2874–2890.

(66) Fey, G.; Reiss, M.; Kersten, H. Interaction of tetracyclines with ribosomal subunits from Escherichia coli. Fluorometric investigation. *Biochemistry* **1973**, *12*, 1160–1164.

(67) ChemAxon Chemicalize. https://chemicalize.com/, 2025; Accessed: 2025-07-30.

(68) Zhang, H.; Shao, X.; Dehez, F.; Cai, W.; Chipot, C. Modulation of Membrane Permeability by Carbon Dioxide. *Journal of Computational Chemistry* **2020**, *41*, 421–426.
(69) Smith, D. J.; Klauda, J. B.; Sodt, A. J. Simulation Best Practices for Lipid Membranes [Article v1.0]. *Living Journal of Computational Molecular Science* **2019**, 1.

(70) Pinsuwan, S.; Alvarez-Núñez, F. A.; Tabibi, E. S.; Yalkowsky, S. H. Degradation kinetics of 4-dedimethylamino sancycline, a new anti-tumor agent, in aqueous solutions. *International journal of pharmaceutics* **1999**, *181*, 31–40.

(71) Pinsuwan, S.; Alvarez-Núñez, F. A.; Tabibi, S. E.; Yalkowsky, S. Spectrophotometric determination of acidity constants of 4-dedimethylamino sancycline (Col-3), a new antitumor drug. *Journal of pharmaceutical sciences* **1999**, *88*, 535–537.

(72) Chen, H.; Ogden, D.; Pant, S.; Cai, W.; Tajkhorshid, E.; Moradi, M.; Roux, B.; Chipot, C. A Companion Guide to the String Method with Swarms of Trajectories: Characterization, Performance, and Pitfalls. *Journal of Chemical Theory and Computation* **2022**, *18*, 1406–1422.
## Supporting Information Available

<table><caption>Table S1: Macro pK<sub>a</sub> values of TC reported by three previous studies.</caption><thead><tr><th>Reference</th><th>pK<sub>a1</sub></th><th>pK<sub>a2</sub></th><th>pK<sub>a3</sub></th></tr></thead><tbody><tr><td>Stephens, et al.<sup>31</sup></td><td>3.30</td><td>7.68</td><td>9.69</td></tr><tr><td>Leeson et al.<sup>64</sup></td><td>3.33</td><td>7.75</td><td>9.61</td></tr><tr><td>Qiang and Adams<sup>65</sup></td><td>3.32</td><td>7.78</td><td>9.58</td></tr></tbody></table>
<table><caption>Table S2: The effective permeability ($P_{eff}$) of TC and the contribution from the neutral tau-tomer ($\frac{P_N f_N}{P_{eff}}$) as the specific permeability ($P_{s}^{other}$) of species other than TC$_Z$ and TC$_N$ (all microstates carrying a net charge and the neutral TC$_Z$, tautomer) changes. The calculation was performed at pH = 6 using the pH-partitioning scheme.</caption><thead><tr><th>P<sub>s</sub> (others)</th><th>P<sub>eff</sub> (cm/s)</th><th>P<sub>N</sub>f<sub>N</sub>/P<sub>eff</sub> (%)</th></tr></thead><tbody><tr><td>0</td><td>3.17 × 10<sup>-5</sup></td><td>99.95%</td></tr><tr><td>P<sub>Z</sub></td><td>3.17 × 10<sup>-5</sup></td><td>99.95%</td></tr><tr><td>10P<sub>Z</sub></td><td>3.17 × 10<sup>-5</sup></td><td>99.93%</td></tr><tr><td>100P<sub>Z</sub></td><td>3.17 × 10<sup>-5</sup></td><td>99.78%</td></tr><tr><td>1000P<sub>Z</sub></td><td>3.22 × 10<sup>-5</sup></td><td>98.24%</td></tr><tr><td>10000P<sub>Z</sub></td><td>3.72 × 10<sup>-5</sup></td><td>85.11%</td></tr><tr><td>100000P<sub>Z</sub></td><td>8.69 × 10<sup>-5</sup></td><td>36.44%</td></tr></tbody></table>
<table><caption>Table S3: The effective permeability of TC ($P_{eff}$) at various pH values obtained using the pH-partitioning scheme (unit: cm/s).</caption><thead><tr><th>pH</th><th>P<sub>eff</sub> (cm/s)</th><th>P<sub>N</sub>f<sub>N</sub> / P<sub>eff</sub> (%)</th></tr></thead><tbody><tr><td>0</td><td>3.14 × 10<sup>-8</sup></td><td>48.0%</td></tr><tr><td>1</td><td>1.67 × 10<sup>-7</sup></td><td>90.2%</td></tr><tr><td>2</td><td>1.50 × 10<sup>-6</sup></td><td>98.88%</td></tr><tr><td>3</td><td>1.03 × 10<sup>-5</sup></td><td>99.84%</td></tr><tr><td>4</td><td>2.66 × 10<sup>-5</sup></td><td>99.94%</td></tr><tr><td>5</td><td>3.16 × 10<sup>-5</sup></td><td>99.95%</td></tr><tr><td>6</td><td>3.17 × 10<sup>-5</sup></td><td>99.95%</td></tr><tr><td>7</td><td>2.74 × 10<sup>-5</sup></td><td>99.94%</td></tr><tr><td>8</td><td>1.14 × 10<sup>-5</sup></td><td>99.86%</td></tr><tr><td>9</td><td>1.41 × 10<sup>-6</sup></td><td>98.84%</td></tr><tr><td>10</td><td>6.88 × 10<sup>-8</sup></td><td>76.27%</td></tr><tr><td>11</td><td>1.70 × 10<sup>-8</sup></td><td>4.17%</td></tr></tbody></table>
<table><caption>Table S4: Physicochemical properties of permeants investigated in this work as given by Chemicalize.⁶⁷</caption><thead><tr><th></th><th>Methanol<br>CH<sub>3</sub>OH</th><th>Urea<br>CH<sub>4</sub>N<sub>2</sub>O</th><th>DDA<br>C<sub>10</sub>H<sub>13</sub>N<sub>5</sub>O<sub>2</sub></th><th>Tetracycline<br>C<sub>22</sub>H<sub>24</sub>N<sub>2</sub>O<sub>8</sub></th></tr></thead><tbody><tr><td>MW</td><td>32</td><td>60</td><td>235</td><td>444.4</td></tr><tr><td>nHDon</td><td>1</td><td>2</td><td>2</td><td>6</td></tr><tr><td>nHAcc</td><td>1</td><td>1</td><td>6</td><td>9</td></tr><tr><td>nC</td><td>1</td><td>1</td><td>10</td><td>22</td></tr><tr><td>nAR</td><td>0</td><td>0</td><td>2</td><td>4</td></tr><tr><td>ClogP</td><td>-0.52</td><td>-1.36</td><td>-0.48</td><td>-3.50</td></tr></tbody></table>
<table><caption>Table S5: Simulation times of WTM-eABF calculations, the free energy barriers and wells of four permeants across POPC bilayers with different sizes.</caption><thead><tr><th>Permeant</th><th>Bilayer</th><th>Time (μs)</th><th>w<sub>peak</sub>(z) (kcal/mol)</th><th>w<sub>valley</sub>(z) (kcal/mol)</th><th>Valley Position (Å)</th></tr></thead><tbody><tr><td rowspan="4">Methanol</td><td>32</td><td>3.0</td><td>3.4 ± 0.0</td><td>-0.0 ± 0.0</td><td>± 28.5</td></tr><tr><td>64</td><td>3.0</td><td>3.6 ± 0.2</td><td>-0.0 ± 0.2</td><td>± 30.9</td></tr><tr><td>128</td><td>4.0</td><td>3.6 ± 0.1</td><td>-0.0 ± 0.1</td><td>± 37.7</td></tr><tr><td>256</td><td>3.5</td><td>3.6 ± 0.2</td><td>-0.0 ± 0.2</td><td>± 32.3</td></tr><tr><td rowspan="4">Urea</td><td>32</td><td>2.9</td><td>10.4 ± 0.0</td><td>-0.2 ± 0.1</td><td>± 24.7</td></tr><tr><td>64</td><td>3.0</td><td>10.6 ± 0.2</td><td>-0.1 ± 0.2</td><td>± 24.5</td></tr><tr><td>128</td><td>3.0</td><td>10.1 ± 0.1</td><td>-0.1 ± 0.1</td><td>± 25.5</td></tr><tr><td>256</td><td>2.3</td><td>10.0 ± 0.1</td><td>-0.1 ± 0.1</td><td>± 24.5</td></tr><tr><td rowspan="4">DDA</td><td>32</td><td>3.0</td><td>7.7 ± 0.1</td><td>-0.1 ± 0.1</td><td>± 29.1</td></tr><tr><td>64</td><td>6.0</td><td>7.7 ± 0.1</td><td>-0.1 ± 0.1</td><td>± 28.7</td></tr><tr><td>128</td><td>3.0</td><td>8.1 ± 0.0</td><td>-0.1 ± 0.0</td><td>± 31.3</td></tr><tr><td>256</td><td>3.3</td><td>7.4 ± 0.1</td><td>-0.0 ± 0.1</td><td>± 35.5</td></tr><tr><td rowspan="4">TC<sub>N</sub></td><td>32</td><td>5.0</td><td>8.1 ± 0.5</td><td>-4.9 ± 0.4</td><td>± 14.3</td></tr><tr><td>64</td><td>9.0</td><td>7.0 ± 0.6</td><td>-5.0 ± 0.6</td><td>± 13.7</td></tr><tr><td>128</td><td>8.0</td><td>4.4 ± 0.1</td><td>-4.9 ± 0.2</td><td>± 14.1</td></tr><tr><td>256</td><td>6.8</td><td>2.2 ± 0.7</td><td>-4.3 ± 0.7</td><td>± 13.5</td></tr><tr><td>TC<sub>Z</sub></td><td>128</td><td>9.0</td><td>13.6 ± 0.2</td><td>-2.1 ± 0.5</td><td>± 18.1</td></tr></tbody></table>
<table><caption>Table S6: Effective permeability ($P_{eff}$) of TC according to the pH-partitioning scheme as a different model compound, 4-Dedimethylamino sancycline (Col-3), is used to estimate $pK_{a1}$. The macroscopic $pK_{a1}$ values of Col-3 given by Chemicalize<sup>67</sup> and Pinsuwan et al.<sup>70,71</sup> are used.</caption><thead><tr><td>pK<sub>21</sub></td><td>P<sub>eff</sub> (cm/s)</td><td>log P<sub>eff</sub></td></tr></thead><tbody><tr><td>5.82<sup>67</sup></td><td>2.25 × 10<sup>-5</sup></td><td>-4.65</td></tr><tr><td>5.64<sup>70</sup></td><td>1.50 × 10<sup>-5</sup></td><td>-4.83</td></tr><tr><td>5.40<sup>71</sup></td><td>8.64 × 10<sup>-6</sup></td><td>-5.06</td></tr></tbody></table>
<table><caption>Table S7: Effective permeability ($P_{eff}$) of TC according to the pH-partitioning scheme obtained by assigning the macroscopic $pK_{a2}$ and $pK_{a3}$ to the microscopic $pk_2$ and $pk_3$, respectively (or in the reverse order). The macroscopic $pK_{a1} = 3.33$ is assigned to $pk_1$ in both cases. The fractions of $TC_N$, $TC_Z$, and $TC_Z'$ are estimated using thermodynamic cycles involving $A^0B^+C^0$, $A^-B^+C^0$, $A^0B^0C^0$ and $A^0B^+C^-$ in Fig. 2.</caption><thead><tr><td>Assignment</td><td>P<sub>eff</sub> (cm/s)</td><td>log P<sub>eff</sub></td></tr></thead><tbody><tr><td>pk<sub>2</sub> = 7.75, pk<sub>3</sub> = 9.61</td><td>8.43 × 10<sup>-7</sup></td><td>-6.07</td></tr><tr><td>pk<sub>2</sub> = 9.61, pk<sub>3</sub> = 7.75</td><td>2.79 × 10<sup>-8</sup></td><td>-7.55</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>Figure S1: Convergence of the PMF calculation of TC<sub>Z</sub> and TC<sub>N</sub> in a 128-POPC bilayer assessed by the RMSD of the potential gradients ∇w(z) and the evolution of the PMF w(z) as the simulation time increases.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure S2: (A) The fractional diffusivity profiles of TC<sub>Z</sub> and TC<sub>N</sub> in a 128-POPC bilayer. (B) Profiles of the estimated fractional order α.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure S3: Lipid phosphorus atom occupancy at an isovalue of 0.02 for TC~N~ in POPC bilayers of different sizes, calculated from trajectories with TC~N~ in the region $-5 < z < 5$ Å.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure S4: Convergence of PMF profiles for TC<sub>N</sub> (A), methanol (B), urea (C), and DDA (D) in POPC bilayers of varying sizes.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure S5: RMSD of the potential gradients ($\nabla w(z)$) as a function of simulation time $t$ for the four permeants in membrane systems containing 32, 64, 128, and 256 POPC, respectively.</figcaption></figure>
<figure><img src="image_20.png" /></figure>

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

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

<figure><img src="image_23.png" /><figcaption>Figure S6: Membrane thickness as a function of the lateral distance in the *xy* plane from the center of each permeant (*r*) in POPC bilayers of different sizes, shown for (A) TC<sub>N</sub>, (B) methanol, (C) urea, and (D) DDA.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure S7: The membrane thickness heatmaps of TC<sub>N</sub>, methanol, urea, and DDA when the permeants are located within the region -5 < z < 5 in POPC bilayers of different sizes.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Figure S8: (A) Orientation of TC<sub>N</sub>, as indicated by the angle between its vector $\vec{v}_2$ and the z-axis, with error bars representing standard deviations. (B) Probability density distributions of the orientation of TC<sub>N</sub> in the regions $-1 \le z \le 1$ Å and $3 \le z \le 5$ Å of POPC bilayers with different sizes.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure S9: Probability density maps of the angle between $\vec{v}_2$ and the z-axis of TC$_N$ in bilayers containing 32, 64, 128, and 256 POPC, respectively.</figcaption></figure>