Injecting Carbon Nanostructures in Living Cells Anastasios Gotzias∗ a.gotzias@inn.demokritos.gr Institute of Nanoscience and Nanotechnology Athens, Greece ABSTRACT most well studied approach to delve into such processes is to use Carbon nanoparticles are currently proposed as reinforcing agents molecular dynamics. However, lengthy simulation runs of “brute- in synthetic biological membranes, able to be embedded in liv- force” molecular dynamics, typically on the nanosecond time scale, ing cells and membrane bilayers. Within a biological environment, would be inefficient to capture the long-time scales of typical bio- porous carbons are anticipated to carry out specific actions, similar logical events, which are frequently on the microsecond or millisec- to the functionality of known assemblies of biological channels like ond time scale. More important, the dissociation of nanoparticles cyclic peptides and aquaporins. An attainable approach to delve through interfaces of cosolvents and bilayers obtains high free en- into the mechanism of how carbon pass through the lipid matrices ergy barriers that cannot be explored using conventional sampling is to use molecular dynamics (MD) simulations. The mechanism methods.[18] This is because, the probability that a spontaneous consists of different stages, the relative free energies of which may fluctuation will bring the system on top of the barrier would be lie far apart in phase space. This induces high energy barriers be- vanishingly small.[7, 16] tween the stages, that cannot be crossed in a single simulation. Such These challenges can be addressed through application of spe- observations are addressed through the application of multi-stage cialized sampling techniques such as umbrella sampling and adap- workflows, where we utilize explicit sampling schemes in every tive force biasing. Such techniques usually require a predefined stage, ranging form grand canonical partitions, for the loading number of executions of single computational tasks. A series of and release of drug substances, to pulling and umbrella sampling advanced sampling techniques can be algorithmically combined in simulations, for the dissociation of nanoparticles. The successful de- multi - stage workflows, to handle complex and highly demand- velopment of workflows relies on the encoding of the dependencies ing computational processes, like those enrolled in bio molecular between the stages and the tasks and the assurance that data and simulations.[5, 19] Arguably, nowhere is the importance of work- parameter variables move between the multi - stage components, flows greater than in biomolecular sciences where the scientific appropriately. The scope is to use the workflow as a descriptor to outcome is intricately intertwined with the ability to execute work- train machine learning models for parameter verification and free flows and computational campaigns successfully.[1] energy calulation methods for carbon - lipid interfaces. 2 METHODOLOGY CCS CONCEPTS We formulate a multi-stage workflow application to encode the • Computing methodologies → Simulation support systems; entire process in which carbon nanoparticles land on, bind to and Molecular simulation. translocate through a lipid environment and release a cargo. This can be accomplished in four sequential computing stages. The first KEYWORDS stage describes the adsorption simulations of the drug substance molecular dynamics, lipid bilayer, porous carbons into the pores of the nanoparticle. The second stage performs the pulling of the nanoparticle into the bilayer (figure 1). The third stage 1 INTRODUCTION is where the nanoparticle is embedded into the membrane and the solvation free energies are computed by decoupling the interfacial If porous carbons are to be exploited as drug delivery systems, it interactions. The forth stage prescribes a model for the diffusivity, is of both fundamental and practical interest to understand how where the cargo substance exits the space of confinement and dis- the carbon interface links to the cholesterol supporters of living sociates to an arbitrary far distance from the nanoparticle. Herefter, cells.[3] Carbons may have nanopores of a size comparable to that we name the different stages of the workflow as 𝑖) adsorption stage, of endogenous protein channels but mimicking their affinity and 𝑖𝑖) pulling, 𝑖𝑖𝑖) decoupling and 𝑖𝑣) drug release stage, respectively. transport properties remains challenging.[6] For instance, surface The four stages of the workflow are partitioned in several sub- functional groups may have adverse effects on the integrity of the tasks (jobs), the development of which, takes place in separate lipid bilayer as they can be toxic.[4, 12, 17] actions. The four stages are then merged in the workflow, i.e., a uni- The entire mechanism of carbon nanoparticles entering into fied module capable to be executed in a single submission. We use and exiting from the lipid environment awaits consensus.[15] The the term "workflow", to express a front-end application handling a ∗ Corresponding author four-stage simulation problem robustly, branching decisions during the stages without the need of user interaction. This development entails the encoding of dependencies of the tasks and stages and AINST2020, September 02–04, 2020, Athens, Greece the assurance that the data and parameter variables move between Copyright © 2020 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). the components and tasks, appropriately. This is important because most of the tasks in a workflow use dependencies from a different Figure 2: Indicative pathways for nanoparticle injection into a membrane bilayer. We create a single pathway (𝛾, 𝛿) us- ing pulling simulation and we perform umbrella sampling along the pathway to explore numerous possible trajectories (Γ𝑖 , Δ𝑖 ) for the injection. The trajectory inside the bilayer (𝛿) is sampled repeatedly, in stage 𝑖𝑖 using umbrella simulations Figure 1: Injecting a single graphene sheet (red hexagonal and in stage 𝑖𝑖𝑖 by decoupling the carbon - lipid interactions, lattice) into a lipid bilayer. The lipids are shown with gray until the estimates for the relative free energy between the lines, the hydrophilic head groups of the lipids are shown two stages statistically converge. with red and blue nodes. Water molecules are shown with cyan points, ions in water are shown with blue spheres. on the drug, that the simulation should definitely take into account. stage and they can be only executed once all of their dependencies Within the development of stage 𝑖𝑣 in the workflow (i.e., drug re- have been completed. Although there have been significant advan- lease stage), we employ a revised caging verification algorithm tages in the state-of-the-theory and practice in workflows, the state that is able to chemically evaluate and predict the hypothetical of workflow development, execution and extension leaves much formation of lipid - substance (host - guest) molecular complexes. scope for improvement. 4 EMBEDDED NANOPARTICLES 3 POROUS CARBONS Many studies that use molecular simulation to describe the pene- From both chemical and technical perspectives, porous carbons tration of membrane cells by carbon sorbents, report contradictory have an important feature; internal cavities. Like in other types of results. Some of these studies depict the lipids attached on the framework materials possessing cavities, substances under confin- carbon surface forming monolayer around the nanoparticle and ment are involved in supramolecular interactions, in particular of blocking the pore channels. Different studies report that the lipids the host-guest type. To discover, whether one substance can access are selective to a particular size of nanoparticles, provided that a specific cavity is a challenging task, because the size and shape their body is hydrophobic.[13] One side of the nanoparticle has of the substance can be very complex. With the nanocarbon model to be shorter than the thickness of the membrane, otherwise the in hand, the only missing component of the theoretical caging pre- nanoparticle leans in a sideways orientation, in order to maximize diction is an algorithm that takes the two geometries as input and its interface contact with the lipids. On the other hand, oxygen determines whether the cavity can encapsulate a substance of an containing functional groups at the rim of the pore channels inter- arbitrary shape. Algorithms of such type are extensivelly used in act with the hydrophilic head groups of the lipid bilayer forming the pore size analysis of crystaline porous solids (metal organic and energetically favorable adsoption sites. The functional moieties zeolitic imidazolate frameworks), where these solids are evaluated on the carbon surface, especially the highly polar ones are great as selective gas filters.[8–10] However, compared to zeolite - type contributors of the insertion process. The polar groups affect the solids, membrane bilayer simulations can depict different, more potential mean force of the membrane penetration so radically that intricate caging complexes. The time the cargo substance escapes a they may render the membrane impermeable to the nanoparticle. pore channel, it can be encapsulated by the lipid macromolecules. However, in most simulation studies, carbons are initilally embed- The lipids configure a cage-like cavity around the drug, that appears ded in the lipid bilayer without any description of how they have like a molecular trap. The trap imposes strong position restraints reached that place. Most studies also employ a unified force field, although it is argued that the surface functional groups should be REFERENCES interpreted with explicit interaction formulas. [2, 20] [1] Kira A. Armacost, Sereina Riniker, and Zoe Cournia. 2020. Novel Direc- We set a reversible path in the 𝑃-𝑇 plane that connects the cur- tions in Free Energy Methods and Applications. Journal of Chemical Infor- mation and Modeling 60, 1 (2020), 1–5. https://doi.org/10.1021/acs.jcim.9b01174 rent simulation system with some reference system of known free arXiv:https://doi.org/10.1021/acs.jcim.9b01174 PMID: 31983210. energy. This prescription implies the use of pulling simulations [2] Mert Atilhan, Luciano T. Costa, and Santiago Aparicio. 2019. On the interaction between carbon nanomaterials and lipid biomembranes. Journal of Molecular (stage 𝑖𝑖, in the workflow).[11] In pull codes, we apply a constant Liquids 295 (2019), 111714. https://doi.org/10.1016/j.molliq.2019.111714 force spring along a reaction coordinate (path), to gradually displace [3] Monica Boffito, Rossella Laurano, Dimitra Giasafaki, Theodore Steriotis, Athana- the nanoparticle from a reference point (point A) to an arbitrary sios Papadopoulos, Chiara Tonda-Turo, Claudio Cassino, Georgia Charalam- bopoulou, and Gianluca Ciardelli. 2020. Embedding Ordered Mesoporous location inside the bilayer (point B), that is the system of interest Carbons into Thermosensitive Hydrogels: A Cutting-Edge Strategy to Vehic- (figure 2). We compute the derivatives of the free energy on the con- ulate a Cargo and Control Its Release Profile. Nanomaterials 10, 11 (2020). secutive steps of this path and integrate. The system of interest may https://doi.org/10.3390/nano10112165 [4] Junlang Chen, Guoquan Zhou, Liang Chen, Yu Wang, Xiaogang Wang, and differ from the reference system, not only in its thermodynamic Songwei Zeng. 2016. Interaction of Graphene and its Oxide with Lipid Mem- state variables but also in its Hamiltonian. This makes possible a brane: A Molecular Dynamics Simulation Study. The Journal of Physical Chem- istry C 120, 11 (2016), 6225–6231. https://doi.org/10.1021/acs.jpcc.5b10635 much wide variety of reference systems and reversible paths. This arXiv:https://doi.org/10.1021/acs.jpcc.5b10635 approach is followed in the stage 𝑖𝑖𝑖 of the workflow where we [5] Clara D. Christ and Wilfred F. van Gunsteren. 2007. Enveloping distribution make an alchemical change on the system of interest.[14] sampling: A method to calculate free energy differences from a single simulation. The Journal of Chemical Physics 126, 18 (2007), 184110. https://doi.org/10.1063/1. In the decoupling stage (stage 𝑖𝑖𝑖), we use point B from the pulling 2730508 arXiv:https://doi.org/10.1063/1.2730508 stage as the initial configuration. We remove the nanopartice from [6] Mariangela Fedel. 2020. Hemocompatibility of Carbon Nanostructures. C — the solvent by varying a decoupling parameter 𝜆 ∈ [0, 1] in steps Journal of Carbon Research 6, 1 (Mar 2020), 12. https://doi.org/10.3390/c6010012 [7] Daan Frenkel and Berend Smit. 1996. Understanding molecular simulation : from of 𝑑𝜆. The decoupling parameter 𝜆, parameterizes the atomistic algorithms to applications. 2nd ed. Vol. 50. https://doi.org/10.1063/1.881812 interactions between the solvent and the nanoparticle. 𝜆 = 0, cor- [8] Anastasios Gotzias. 2017. The effect of gme topology on multicomponent ad- sorption in zeolitic imidazolate frameworks. Phys. Chem. Chem. Phys. 19 (2017), responds to the state where the interactions are full (point B) and 871–877. Issue 1. https://doi.org/10.1039/C6CP06036F 𝜆 = 1, corresponds to the state where the nanoparticle does not [9] Anastasios Gotzias. 2019. Calculating adsorption isotherms using Lennard Jones interact with the lipids as if it is simulated in vacuum (point C).[16] particle density distributions. Journal of Physics: Condensed Matter 31, 43 (jul 2019), 435901. https://doi.org/10.1088/1361-648x/ab2c94 This involves the execusion of independent molecular dynamics [10] Anastasios Gotzias, Evangelos Kouvelos, and Andreas Sapalidis. 2018. Computing simulations for the different values of 𝜆. From the simulations we the temperature dependence of adsorption selectivity in porous solids. Surface get the average derivative of the parameterised Hamiltonian. Then and Coatings Technology 350 (2018), 95 – 100. https://doi.org/10.1016/j.surfcoat. 2018.07.003 we compute the free energy, Δ𝐺 𝐵𝐶 , using integration. [11] Anastasios Gotzias and Andreas Sapalidis. 2020. Pulling Simulations and Hy- The pulling and decoupling stages use the same system of interest drogen Sorption Modelling on Carbon Nanotube Bundles. C 6, 1 (2020), 11. https://doi.org/10.3390/c6010011 (point B), while their reference systems differ only on the type of [12] Ruohai Guo, Jian Mao, and Li-Tang Yan. 2013. Computer simulation of cell solvent, that is water in point A for the pulling stage and vacuum entry of graphene nanosheet. Biomaterials 34, 17 (2013), 4296 – 4301. https: in point C for the decoupling stage. We can compute the difference //doi.org/10.1016/j.biomaterials.2013.02.047 [13] Georgia Karataraki, Andreas Sapalidis, Elena Tocci, and Anastasios Gotzias. on the free energy change between points A and C by summing 2019. Molecular Dynamics of Water Embedded Carbon Nanocones: Surface the free energy changes of the two stages, Δ𝐺𝐴𝐶 = Δ𝐺𝐴𝐵 + Δ𝐺 𝐵𝐶 . Waves Observation. Computation 7, 3 (2019), 50. https://doi.org/10.3390/ The next step is to make a subtle change on an input parameter (i.e., computation7030050 [14] David L Mobley, John D Chodera, and Ken A Dill. 2006. On the use of orientational the value of a Lennard jones parameter) and run the pulling and restraints and symmetry corrections in alchemical free energy calculations. The decoupling stages. We change again this value and run this process Journal of chemical physics 125, 8 (08 2006), 084902–084902. https://doi.org/10. 1063/1.2221683 iteratively, until the free energy change, Δ𝐺𝐴𝐶 , converges. [15] Elle Puigpelat, Jordi Ignés-Mullol, Francesc Sagués, and Ramon Reigada. 2019. Interaction of Graphene Nanoparticles and Lipid Membranes Dis- playing Different Liquid Orderings: A Molecular Dynamics Study. Lang- muir 35, 50 (2019), 16661–16668. https://doi.org/10.1021/acs.langmuir.9b03008 arXiv:https://doi.org/10.1021/acs.langmuir.9b03008 PMID: 31750663. [16] Michael R. Shirts, Jed W. Pitera, William C. Swope, and Vijay S. Pande. 2003. 5 CONCLUDING REMARKS Extremely precise free energy calculations of amino acid side chain analogs: Comparison of common molecular mechanics force fields for proteins. The The role of molecular simulation studies are to discover the key Journal of Chemical Physics 119, 11 (2003), 5740–5761. https://doi.org/10.1063/1. factor at the nanoscale which is usually ignored and provide an 1587119 arXiv:https://doi.org/10.1063/1.1587119 [17] Martin Vögele, Jürgen Köfinger, and Gerhard Hummer. 2018. Molecular dynamics understanding that will break the conventional way of nanoporous simulations of carbon nanotube porins in lipid bilayers. Faraday Discuss. 209 material design and application. Membrane bilayers are compli- (2018), 341–358. Issue 0. https://doi.org/10.1039/C8FD00011E cated molecular systems with several degrees of freedom and corre- [18] Di Wu. 2010. Understanding free-energy perturbation calculations through a model of harmonic oscillators: Theory and implications to improve lated torsional terms. In order to sufficiently sample such systems the sampling efficiency by molecular simulation. The Journal of Chem- it requires increased computational power and smarter sampling ical Physics 133, 24 (2010), 244116. https://doi.org/10.1063/1.3511703 schemes. Using machine learning, we show that the pathway to arXiv:https://doi.org/10.1063/1.3511703 [19] Di Wu and David A. Kofke. 2005. Phase-space overlap measures. II. Design and accurate and reliable methods to compute the free energies of such implementation of staging methods for free-energy calculations. The Journal systems may be clearer than previously thought. This is especially of Chemical Physics 123, 8 (2005), 084109. https://doi.org/10.1063/1.2011391 arXiv:https://doi.org/10.1063/1.2011391 true in the light of new distributed computing techniques, which [20] Songwei Zeng, Yu Ji, Yue Shen, Ruiyao Zhu, Xiaogang Wang, Liang Chen, and provide the greatly increased computational power needed for both Junlang Chen. 2020. Molecular dynamics simulations of loading and unloading the development of improved parameter sets and the sufficient of drug molecule bortezomib on graphene nanosheets. RSC Adv. 10 (2020), 8744– 8750. Issue 15. https://doi.org/10.1039/D0RA00261E sampling of extended ensemble methods.