The 2020 pandemic has brought the chasm between the rapid spread of infectious diseases and the relative slow pace of scientific discovery to stark reality. This particular event is not the first, nor it will be the last. Since the turn of the century three different coronavirus threats have surfaced and emerging infectious diseases will continue to be a threat for the foreseeable future. Current methods for studying betacoronaviruses, responsible for the 2020 pandemic, are technically and time demanding. Viral isolation from field sample is rarely successful and reverse genetics recovery of recombinant virus is labor intensive and expensive. Similarly, the design and testing of new drugs is an extremely lengthy process. To this end, drug repurposing, where existing compounds that have already been tested safe in humans, are redeployed to combat another disease, has emerged as a critical strategy. As an example, Remdesivir, one of the promising COVID-19 drug candidates at the moment, is a broad-spectrum drug that was originally targeted toward a number of viruses, such as SARS-CoV and MERS-CoV, but also studied to counter Ebola virus infection.
Multiple characteristics define the potency of a drug, as well as its side effects. Among them, bioavailability -the fraction of a compound that reaches the intended target -is one of the hardest to estimate reliably and rapidly. Although several factors contribute to in defining bioavailability, passive permeability through different biological membranes is a critical component since independent of the ability of a drug to bind to its target, the inability to reach the intended target often translates into poor in vivo efficacy and may require higher dosages with higher risks of side effects. 1 The viral life cycle encompasses several crucial steps, starting with the attachment of the virus to the cell and finishing with the release of the progeny virions from the host cell. SARS-CoV-2 targets cells through the viral structural spike protein that binds to an enzyme on the host cell membranes. Following receptor binding, the virus uses host cell receptors and endosomes to enter cells. Once inside, viral polyproteins are synthesized that encode for the replicase-transcriptase complex. The virus then synthesizes RNA and structural proteins leading to completion of assembly and release of viral particles. These viral lifecycle steps provide potential targets for drug therapy and most of the drugs currently in clinical trials in-hibit key components of the infection life-cycle, such as viral entry into the host cell, viral replication, and viral RNA synthesis. 2 These targets are localized to particular subcellular compartments, yet current drug design strategies are focused on general bioavailability and rarely address drug delivery to specific tissues or intracellular compartments. Therefore, knowledge of the permeability of antiviral drugs, with reference to both viral and mammalian cell compartments, is vitally important from both a pharmacokinetics and drug design standpoint.
In an effort to accelerate translation of promising drug candidates, as well as eliminate potentially ineffective or even harmful candidate, we present a science-based model to predict the permeation of drug candidates for SARS-CoV-2 into coronavirus membrane and mammalian membranes.
Over the past decades many models have been developed to estimate drugs' permeability. The most commonly used ones are the empirical "rule of five", 3 QSAR (Quantitative Structure-Activity Relationship) 4 and QSPR (Quantitative Structure-Property Relationship) methods. 5 While QSAR/QSPR methods have shown their efficacy in screening the binding of antiviral drugs and targets, the predicted permeability is at best only accurate to an order of magnitude, mostly due to the incorporation of too many descriptors that are not related to A) Protein colored according to secondary structure (purple for alpha helix, red for beta sheet, and gray for random coils). In the lipid bilayer, carbon are shown in gray, nitrogen in red, oxygen in blue, and phosphate in tan; cholesterol is shown uniformly in orange. Hydrogen atoms omitted for clarity. B) The anchoring region (part of the heptad repeated, the transmembrane region, and the C-tails) of the spike protein in the membrane. C) The anchoring region of the spike protein colored by the residue type. Non-polar residues are colored in white, while polar residues are in green. the permeation problem, and the strong dependence on the quality of the training sets employed. 6 On the other end of the accuracy spectrum, several methods have been proposed over the years that leverage a physical description of the passive permeation process, like the inhomogeneous solubility diffusion model and statetransition model. 7 These models, that currently often involve molecular dynamics (MD) simulations, are either computationally too expensive for a rapid application or suffer from limited accuracy when drugs of various sizes and shapes or membranes with different lipid profiles are involved.
To bridge the gap between accuracy and ease of use, we have recently introduced a model to predict the rates of passive permeation of molecules and nanoparticles (up to several nm) into phosphatidylcholine membranes. 8 The Low-Density Area (LDA) model, which can combine both experimental and computational data, has shown accuracy equivalent to fully atomistic MD methods with a much smaller computational cost. Such a result was achieved observing that the permeation process of small compounds is controlled by the formation of low-density areas on the membrane surface, and by factorizing the permeation process into independent contributions of the membrane thermal and pressure fluctuations and characteristics of the molecules (size, shape and lipophilicity).
Leveraging the LDA model, here we report on the permeation kinetics of promising drugs candidates through coronavirus membranes, with the overarching goal of providing inputs to help the design and selection of new or repurposed drugs. Additionally, since the target of many drugs are located inside mammalian cells, we also report on the permeation kinetics of the same drugs through five mammalian membranes. A schematic of the overall approach and data sources used in this study is shown in Fig. 1 . The computations were obtained using a combination of molecular dynamics simulations in conjunction with an augmented LDA model for membranes. Seventy-nine drugs for SARS-CoV-2 were analyzed, that comprise FDA approved drugs, compounds in clinical trials, and preclinical compounds. 2 The computational results are compared with data available in the literature on bioavailability and show the potential of the model to aid in the assessment of drugs' permeability into a variety of organelles.
Coronaviridae is a family of large, enveloped, positive-stranded RNA viruses. The envelope is derived from the host Endoplasmic Reticulum (ER) membrane and consists of a lipid bilayer, in which structural proteins are anchored. The surface glycoprotein spike that forms large protrusions from the virus surface, binds to the host-cell receptor and mediates viral entry. 9, 10 Since surface characteristics of the biological membranes determine the permeability of drugs, we developed a model to describe the anchoring of the spike protein to the lipid bilayer. The spike contains three segments: a large ectodomain, a single-pass transmembrane anchor, and a short intracellular tail. The region outside the viral core forms alpha helices where coiled-coil interactions, through hydrophobic contact, help stabilize the chains. The tails near the C-terminal are mostly polar with charged residues that form random coils inside the viral core (i.e., the cytoplasmic side after membrane fusion with the host cell): interactions between the charged residues and lipid head groups stabilize the end of the transmembrane region of the spike protein. A diagram of a section of the coronavirus envelope composed of lipid bilayer and the anchored spike protein is shown in Fig. 2 .
Based on existing data and estimated secondary structure, 11, 12 we determined that the transmembrane region of the COVID-19 spike protein is formed by the residues in the #1210-1235 range. The charge-charge interactions among residues 1215:TYR, 1219:GLY, and 1223:GLY contribute to hold the three chains together inside the hydrophobic region of the membrane lipids. The density of the spike protein, used in our simulations, reflects the experimental density of about 100 nm 2 to 200 nm 2 per transmembrane proteins (of any kind), estimated from a viral envelop diameter of 80 nm to 120 nm 13 and an approximate number of 800 transmembrane proteins per virion. 14 To reduce the computational cost, after the initial equilibration run used to validate the stability of the protein anchoring region, we considered only part of the protein (Fig. 2B and C) as the rest of the spike protein is far enough (5 nm or more) to affect the formation and life time of low-density areas on the viral membrane surface.
For the lipid bilayer, slightly higher concentrations of sphingomyelin and phosphatidylinositol were found in the viral envelope, compared to the ER membrane from where coronaviruses acquire their membrane envelope. 15 Drugs' Membrane Permeability The LDA model describes the permeation as a gated entry process, regulated by the chance of the molecule to find a low-density area on the membrane surface and by the solubility of the drug in the hydrophobic phase. The four parameters used in the model separate the contributions of the characteristics of the molecules (size, shape, solubility) from the properties of the membranes (lipid density distribution).
In addition to the viral envelope described before (with and without spike protein), we studied the surface properties of five membranes corresponding to various mammalian compartments (lysosome, plasma, Golgi, mitochondrial and endoplasmic reticulum membranes) using all-atom molecular dynamics simulations. Fig. 3 reports the distribution of the permeation times of the 79 drugs in all the membranes relative to the ER membrane, which was chosen as reference since the permeation in this membrane is the fastest for all drugs. The results show that the inclusion of the transmembrane protein in the COVID-19 membrane has a large impact on drug's membrane permeability. The presence of the transmembrane protein affects the frequency, shape, and lifetime of lowdensity areas in the LDA model. Indeed, when comparing the distribution of permeation times for ER and coronavirus membrane, with and without the spike protein, we observe a ∼40% increase in the permeation time due to the presence of proteins and ∼10% due to the change in lipid composition. This effect is mostly due to the reduced mobility of the lipids near a rigid object similar to what has been observed for C 60 fullerenes in lipid bilayers. 16 While the presence of the spike protein in the COVID-19 membrane results in an increased stiffness, the viral membrane is still generally two times more permeable than plasma membrane and lysosomes (before accessible surface and concentration effects are considered). Although currently very few drugs are required to permeate through viral membranes, this difference can potentially be exploited by other classes of compounds, like nanoparticles, which can amplify the effect of membrane stiffness due to their relative large size. The individual times of permeation of drug candidates grouped by viral baits or, when that is unknown, by biological processes, are shown in Fig. 4 . The results show a variety of timescales going from nanoseconds to microseconds depending on the drug and membrane.
Albeit the bioavailability of a drug is affected by many factors, like metabolism rate and plasma protein binding, 17 we find a negative correlation between the bioavailability of drugs and the permeation time predicted by the LDA model. For example, for drugs within the same group, Lisinopril (Fig. 4B , Cell Entry) is the least permeable compound as well as the compound with the lowest bioavailability (25% based on experimental data 18 ). On the other hand, Captopril, Nafamostat, and Camostat, for which we predict permeation times that are one order of magnitude shorter than Lisinopril, have an estimated bioavailability of 100% 1 . Similar observations can be made for Rapamycin (Fig. 4A, N) , which has lower bioavailability than Sapanisertib and Silmitasertib. While undeniably other factors affect the ability of a drug to reach the intended target, these results show how useful membrane permeability predictions can be.
Interestingly, the permeation times of drugs that are currently identified as most promising candidates against COVID-19 during their clinical trial 19 are among the shortest compared to other drugs with similar targets. For example, drugs related to the disruption of viral translation such as Ternatin-4 20 (Fig. 4B Translation) , Zotatifin 21 (Fig. 4A, Nsp2) and Plitidepsin 22 (Fig. 4B, Translation) , and sigma receptors such as Haloperidol 23 and Compound PB28 24 (Fig. 4A Nsp6 and Orf9c) , Melperone, Cloperastine, 25 and Clemastine (Fig. 4B , Protein Processing), and the female hormone progesterone (Fig. 4B, Viral transcription) are among the most permeable compounds in their respective categories. While a direct comparison among these drugs is not straightforward, as they need to cross different membranes in order to reach their intended target, the time of permeation in their least permeable membrane (i.e., lysosome or plasma) is 40 ns on average with an upper bound of 62 ns, while the remaining drugs require an average of 242 ns and a maximum of 4.6 µs.
In this paper, we present a model to compute the permeability of drugs that target various steps in the lifecycle of infection of SARS-CoV-2. The model is based on the LDA theory developed by our group that describes the permeation process as regulated by low-density area on the membrane surface and solubility of drugs in the hydrophobic phase of the membrane itself. Starting from the virus structure, we developed an atomistic model of the spike protein anchored to the lipid bilayer to study the viral envelope. Molecular dynamics simu-lations of the viral and five mammalian membranes were carried out to determine the dynamics and characteristics of the membranes. The computational results were compared with data present in the literature on bioavailability, and show that the LDA model can be applied to predict the permeability of drugs through various membranes and to gain insights on the differences between membranes and potentially organs or organisms. From a drug design, as well as health assessment perspective, this predictive capability can facilitate compound evaluation by ruling out candidates with low permeability or excessive permeability in unintended organs or organelles. Indeed, the LDA model can be used to estimate the absorption selectivity of different targets (e.g., carcinogenic cells, organs) by molecules and nanoparticles without need of large set of preexisting data, making it an invaluable tool for the comprehensive and fast validation pipeline of new compounds.
While the results reported in this paper focus on drugs for SARS-CoV-2, the method is general and flexible and can be employed, with other tools to shorten the effective response time to a global pandemic. Our work provides a means for scientific discovery to move closer to the speed of pandemic spreads and it can be generalized to future events.
We used NAMD 26 version 2.13 with CHARMM Force Field 27 version 36 for lipids, and TIP3P 28 for water molecules. All MD simulations were performed with a timestep of 2 fs while keeping all the C -H and O -H bonds rigid via the SHAKE algorithm. 29 Long range electrostatic interactions were modeled with the particle mesh Ewald method 30 using a 0.1 nm grid spacing, a tolerance of 10 −6 and cubic interpolation; cubic periodic boundary conditions were applied. Temperature was kept constant at 310 K using a Langevin thermostat 31 with a time constant of 1 ps, while pressure was kept at 101.325 kPa using Langevin piston method 32, 33 with period of 50 fs and 25 fs decay. To account for the intrinsic anisotropy of the system, the production simulations were performed in the NPsT ensemble, where the x and y dimensions of the periodic system (coplanar with the membrane) are allowed to vary independently from the z dimension. Non-bonded short-range interactions where smoothly switched to 0 between 1 and 1.2 nm with a X-PLOR switch function.
The lipid compositions for cellular organelles (plasma, lysosome, Golgi, mitochondrial) were taken from Horvath et al., 34 and the ER and viral membrane from van Genderen et al. 15 Table S1 provides more details about the composition of each membrane. Membranes without spike proteins had an approximate size of 8 nm by 8 nm, subject to round up to satisfy the required imposed by the lipid composition. The membranes with spike proteins were extended to guarantee at least 4 nm of separation between the closest atom of the spike protein with its nearest periodic image, resulting in membranes of approximately 15.4 nm by 15.4 nm for the system with the entire protein and 9.5 nm by 9.5 nm for the system with the anchoring region only. System were simulated in explicit solvent to guarantee at least 6 nm of separation between the membrane or protein and their periodic image from top and bottom. The lipids, protein, water molecules and ions were initially placed by the CHARMM-GUI membrane builder. 35 Each membrane was equilibrated for at least 30 ns in the NPsT ensemble described above. The membranes were assumed to be at equilibrium once the time average of the area per lipid in a NPsT simulations was varying less than 10% in 10 ns. The surface properties (LDAs) were then obtained by running a additional 100 ns.
The main structure (residues 27 to 1146) of the spike protein can be found in the protein data bank #6VSB, based on the cryo-EM structure of Wrapp et al. 11 The secondary structure of residue 1146 to 1273 was instead based on the prediction of the C-I-TASSER model. 12 The secondary structure of the alpha helices of the lower chunk was left unmodified, while the random coils were twisted to align the alpha helices perpendicularly to the membrane. Simulations of 30 ns show that the anchor is stable with no lipid rearrangement or vertical translation relative to the membrane.
To better factor the drugs' flexibility in the model, their size/shape was estimated by considering a large sample of conformations at room temperature in water. For this purpose, we run molecular dynamics simulations using the CHARMM General Force Field 36 to model the atomic forces and collected the molecule's conformations every 20 ps. CGenFF 37 software was used to assist in the production of the needed topology files. Each drug molecule was presented at the center of an cubic box of 6 nm in each dimension, containing explicit solvent of water and 0.15 M NaCl. MD setup for drug simulation was the same as the membrane simulations, except the use of an isotropic ensemble (NPT). Each system was minimized for 1000 steps, following 0.5 ns of relaxation and 0.5 ns of production. Every 20 ps the molecule shape was estimated using the protocol described previously 8 and then used to calculate the P LDA for the specific conformation. The final value for P LDA was computed by averaging over the values obtained for all the conformations. Finally, the partition coefficients of the drugs molecules in octanol and water were taken from PubChem website using experimental values when available, or prediction from XLogP3 model 38 otherwise. based permeability model for drug-like compounds. Bioorg. Med. Chem. 2011, 19, 2615-2624.