Not so rigid capsids based on cyclodextrin complexes: keys to design

Not so rigid capsids based on cyclodextrin complexes: Keys to design Fabián Suárez-Lestón a , Pablo F. Garrido a , Ángel Piñeiro a, ⇑ , Rebeca Garcia-Fandino b, ⇑ a Departamento de Física Aplicada, Facultade de Física, Universidade de Santiago de Compostela, E-15782 Santiago de Compostela, Spain b Departamento de Química Orgánica, Centro Singular de Investigación en Química Biolóxica e Materiais Moleculares (CiQUS), Universidade de Santiago de Compostela, Campus Vida s/n, E-15782 Santiago de Compostela, Spain graphical abstract article info Article history: Received 30 March 2022 Revised 14 May 2022 Accepted 16 May 2022 Available online 19 May 2022 Keywords: Cyclodextrins Membranes Monolayers Bilayers Molecular dynamics simulations Self-assembly abstract Hypothesis: Membranes based on cyclodextrin complexes can be used as functional nanocarrier envelopers by chemical modifications of the cyclodextrin hydroxyl groups or by encapsulating different ligands in their cavities. Experiments: Molecular dynamics simulations of monolayers and bilayers based on supramolecular complexes consisting of two a or b-cyclodextrin and one sodium dodecylsulfate or dodecane at 283 K and at 298 K were performed. Findings: It is shown that the structure and main interactions stabilizing the membranes, as well as their permeability to water and ions can be tuned by changing the cyclodextrin, the ligand, the number of layers or/and the temperature. These results provide new evidences about both their dynamic nature and the interactions responsible for the stabilization of the membranes and will facilitate the design of new functional capsides and applications based on cyclodextrin complexes. Ó2022 The Authors. Published by Elsevier Inc. Thisis an openaccess article underthe CCBY license(http:// creativecommons.org/licenses/by/4.0/). 1. Introduction Lipid-based membranes are the most typical structures employed by nature to define the boundary between living organisms and their environment. These membranes are formed by two monolayers of lipid molecules, faced to each other, with the lipid https://doi.org/10.1016/j.jcis.2022.05.098 0021-9797/Ó2022 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ⇑ Corresponding authors. E-mail addresses: [email protected] (Á. Piñeiro), rebeca.garcia.fandino@usc. es (R. Garcia-Fandino). Journal of Colloid and Interface Science 623 (2022) 938–946 Contents lists available at ScienceDirect Journal of Colloid and Interface Science journal homepage: www.elsevier.com/locate/jcis polar heads oriented to the aqueous media in both leaflets. Lipid membranes are spontaneously formed when a high enough concentration of these molecules is present in water, driven by the so-called hydrophobic interactions. In general, most biological membranes are formed by a suitable combination of a large variety of lipids, proteins and other heterogeneous molecules that affect their structure as well as their mechanical and dynamical properties. Given the importance and intriguing properties of biological membranes there is a strong motivation to build biomimetic systems using non-living precursors. Artificial membranes have been built based on different synthetic components, including polymers and surfactants [1]. They are typically employed to design new materials or for biotechnological applications [2]. Recently, a new type of membranes based on cyclodextrins (CDs) that selfassemble at polar-nonpolar interfaces has been proposed [3–7]. CDs are well known oligosaccharides formed by 6 ( a -), 7 (b-) or 8( c -) a - D -glucopyranoside units. Yang et al. [8] revealed that CD complexes with 2:1 stoichiometry formed by two b-cyclodextrins and one sodium dodecylsulfate (SDS) molecule (SDS 1 bCD 2 ) spontaneously form giant capsids, tubular or multilamellar structures, depending on the concentration, in aqueous solution. These authors proposed the use of such compartmentalized and rigid structures to design new functional materials able to mimic and to extend the applications of protein capsids. In addition to a detailed structural characterization, they hypothesized the presence of specific interactions to stabilize the membranes, including a network of lateral hydrogen bonds between contiguous SDS 1 - bCD 2 complexes. In contrast, the entropy increase due to the release of a large number of highly ordered water molecules has been recently proposed to explain the adsorption of SDS 1 a CD 2 complexes to polar-nonpolar interfaces [9]. Given the connection between such adsorption process and the self-assembly in aqueous solution, the same type of interactions are expected to drive both phenomena. Thus, these systems are not well understood yet and they require further investigation. In addition to the interactions leading to and stabilizing these systems, their potential abilitiy to filter water, ions or other molecules has not been studied yet. The existence of membranes based on CDs proposed by Yang et al. [8] is specially interesting because, in addition to be employed to encapsulate a variety of molecular systems in the self-assembled compartments spontaneously formed by them, they can be easily funcionalized by introducing different groups in their cavities or by substituting the hydroxyl groups of the CDs. Thus, specific structures with a funcionalized surface that can simultaneously act as molecular carriers, can be designed based on CDs. The understanding of the interactions that govern the formation of such membranes, as well as their ability to permeate different molecules, are key to design new applications. In this work, we focus precisely in these features, as well as in the chemical nature of the building blocks, by a computational approach to these systems. Our study is based on relatively long (2 l s) molecular dynamics (MD) trajectories of membranes formed by SDS 1 - bCD 2 complexes at 283 K and at 298 K. These temperatures were selected because it was shown that the ratio between 1:1 and 2:1 complexes dramatically changes in this interval [10,11]. Simulations of monolayers and bilayers based on 2:1 complexes as building blocks were also performed to assess the possible stability of both types of assemblies as well as their different behavior. Three additional similar systems, based on complexes with 2:1 stoichiometry formed by 2 a CD and 1 SDS or dodecane (DOD), and by 2bCD and 1 DOD, were also studied at the same temperatures (Fig. 1). Both monolayers and bilayers formed by such building blocks were considered. This allows to assess the impact of the charge in the guest molecule as well as to compare the structure and properties of membranes formed by a CD with those formed by bCD and of monolayers vs bilayers. The structure and stability of the possible membranes as well as their ability to permeate water and Na + ions were characterized. 2. Methods 2.1. Set up of the simulation boxes and MD simulation parameters Manually preassembled X 1 a CD2 and X 1 bCD 2 complexes, with X being SDS or DOD, were introduced in 5x5x5 nm 3 cubic boxes solvated with approximately 4000 SPC [12] water molecules. For the systems containing SDS, one sodium counterion was added. The resulting systems were minimized using the steepest descent algorithm, with the 54a7 version of the GROMOS force field [13] and the GROMACS 2018 package [14,15]. Then, 5-ns-long MD trajectories were generated at 298 K and 1 bar to equilibrate the systems. The V-rescale thermostat [16] with a coupling constant of 0.1 ps and an isotropic Parrinello-Rahman barostat [17] with a compressibility of 4.510 5 bar 1 were employed in these simulations. The topology of cyclodextrins and SDS were taken from previous works [9,18–20] and the topology of DOD was adapted from that of SDS. The optimized structures of the simulated 2:1 complexes were rotated to have their symmetry axis parallel to the Z axis, and replicated in the XY plane to generate 5x5 monolayers. For the systems based on a CD hexagonal lattices with a = b = 1.3 8 nm were built while for those with bCD rombic crystals with a = b = 1.52 nm and an angle of 104ºwere generated, as suggested from experimental results [2]. The resulting arrangements were introduced in prism-shaped simulation cells in such a way that the XY surface is adapted to tightly contain the complexes, displayed in a 55 grid with the appropriate lattice. The size of the box in the Z dimension was 10 nm. Additional systems were built by replicating the monolayers in the Z dimension to form bilayers based on the 2:1 complexes. These assemblies were also introduced in prism-shaped boxes with a height of 11 nm. The resulting boxes were solvated with SPC water molecules and Na + counterions were added, when necessary, to neutralize the total charge of each system. The total number of water molecules varied between 11,454 and 16,123 depending on the system. In this way, eight different simulation boxes using a CD or bCD, SDS or DOD, as well as monolayers and bilayers for each combination of molecules, were obtained. After a steepest descent minimization, 2 l s-long molecular dynamics trajectories were obtained for each of these systems at 283 K and at 298 K using a Nose-Hoover thermostat. For these simulations a Parrinello-Rahman semiisotropic pressure control was employed, to allow independent fluctuations of the box dimensions in the XY plane and in the Z dimension. In all MD trajectories, a timestep of 2 fs was employed to integrate the motion equations, using the leap frog algorithm [21]. Long range electrostatic interactions were calculated using the particle mesh Ewald method [22,23] with a real-space cutoff of 1.2 nm, a 0.15 nm spaced grid, and fourth-order B-spline interpolation. The Ewald sum in three dimensions with a correction term (EW3DC) was used to avoid artifacts due to interactions between periodic images in the Z direction. The SETTLE [24] algorithm was employed to constrain the bond lengths and angles of water molecules while the LINCS [25] algorithm was employed to constrain the bond lengths of the CD molecules. During the MD simulations, coordinates and energies were stored every 14 ps for analysis. 2.2. Analysis of the trajectories The area of the simulation boxes was taken as a quantitative indicator to assess the stability of the membranes. Hydrogen bonds between different groups were calculated using the hbonds GROMACS tool: (i) primary-primary, primary-secondary and Fabián Suárez-Lestón, P.F. Garrido, Á. Piñeiro et al. Journal of Colloid and Interface Science 623 (2022) 938–946 939 secondary-secondary hydroxyl groups between CDs (excluding intramolecular interactions); (ii) interactions between monomers of the same 2:1 complex; (iii) interactions between CDs of different monolayers for the bilayer systems; and (iv) interactions between CDs and water molecules. Water density profiles as a function of time, averaged every ns, were computed along the first microsecond and plotted as a greyscale colormap. For the second microsecond of each trajectory, the average density profile of water molecules and oxygen atoms of the CDs joining the glucopyranoside rings (labelled as O4) was also obtained and represented over the same plot. The transport of water and ions across the membranes was also quantified by following the Z coordinate of every single particle along the whole MD trajectories. Diffusion coefficients were determined from the fitting of the mean lateral displacement distributions to the two-dimensional random walk equation [26] over the last l s of the trajectories using time windows of 1 ns and 5 ns. The global arrangement of the molecules has also been independently investigated for each monolayer, even in the systems with two layers of X 1 CD 2 structures. To overcome the perturbations of the atoms lattice and to take into account its time evolution, Voronoi tessellation has been chosen to characterise the topology [27]. This analysis was based on the quantification of the angles consisting the Voronoi cells obtained from the projection of the average position of the O4 atoms of each complex in the XY plane. The equivalence with the square, rhombic and hexagonal lattices, together with the perturbed square, is shown in the Fig. 2. The Voronoi cells were obtained considering the periodic boundary conditions to avoid border artifacts. The angle of all the Voronoi vertices, as well as the the distances between them, were determined as a function of time every 1.4 ns (i.e., every100 frames). This analysis was performed by using a segment of 500 ns taken from each trajectory and all the results were plotted together as a distribution. The fragment of the trajectory used in each case was taken to avoid significant perturbations that eventually appear in some monolayers when some complex rotates, disaligning its symmetry axis from the Z axis. All the analysis was performed using GROMACS commands combined with specifically developed Python scripts using mainly the MDAnalysis [28,29], Numpy [30], Matplotlib [31] and Scipy [32] packages. 3. Results Monolayers and bilayers based on a CD or bCD with SDS or DOD forming complexes with 2:1 stoichiometry (SDS 1 a CD 2 , SDS 1 bCD 2 , DOD 1 a CD 2 , DOD 1 bCD 2 ) were built and simulated for 2 l sat 283 K and at 298 K, giving a total of sixteen MD trajectories. In general, all the planar assemblies were stable although significant perturbations and different structural and dynamic behaviors were observed depending on the temperature of the system, on the molecules consisting the complexes and on the number of layers forming the membrane. Snapshots of the initial and final structures as well as the result of the analysis described in the methods section for all the systems are included in the Supplementary Informa- Fig. 1. A. Schematic representation of the complexes studied in this work: SDS 1 a CD 2 , SDS 1 bCD 2 , DOD 1 a CD 2 , DOD 1 bCD 2 forming monolayers (B) or bilayers (C). CDs that are closer and further from the SDS ionic head are represented in orange and yellow colors, respectively. Fig. 2. Voronoi cells (green) of four different geometric configurations for the selected atoms (white circles). Fabián Suárez-Lestón, P.F. Garrido, Á. Piñeiro et al. Journal of Colloid and Interface Science 623 (2022) 938–946 940 tion (SI) file (Figs S1–S139 and Tables S1–S7). The most important features of the studied systems will be described in what follows: 3.1. Lateral packing of the complexes For the systems with a CD, the average area available per supramolecular complex in the XY plane was quite stable as a function of time -after the first 200 ns- and similar for the different trajectories (Figs. S121–S124, Table S1 and Fig. 3). The differences observed for both temperatures, for monolayers and bilayers, and using SDS or DOD threading the CDs were very small. The minimum value was observed for the monolayer of SDS 1 a CD 2 at 283 K (174.34 ± 0.11 Å 2 ) while the maximum value was observed for the bilayer of DOD 1 a CD 2 at the same temperature (176.92 ± 0. 10 Å 2 ), (Fig. 3 and Table S1). The area per complex for the simulations with bCD was much larger and noisy than for those using a CD (Figs. S121-S124). This was partially expected since the former CD is wider than the latter, but it seems that the packing of the membranes based on bCD is lower than that of those based on a CD. These simulations showed to be also significantly sensitive to the studied variables. In general, the structures are closer to each other in the XY plane for the monolayer systems than when two layers of complexes are present. The area per complex is not sentitive to the ligand encapsulated in the bCD and the systems at 298 K require more area per complex than those at 283 K. An exception to this behavior is observed for the simulation of a bilayer with SDS at 283 K, which has an abnormally large area per complex (238.96 ± 0.13 Å 2 ) compared to that of the equivalent system at 298 K (234.33 ± 0.15 Å 2 ) and almost equivalent to that of the bilayer formed by DOD 1 bCD 2 complexes at 298 K (238.15 ± 0.19 Å 2 ). This large area is arises from the significant tilting of the complexes in these simulations (see images in the SI) but the reason of this tilting is not clear. The average area per complex obtained for the monolayers of the systems with bCD was very similar for both ligands and clearly temperature dependent: (224.17 ± 0.21 Å 2 and 224.80 ± 0.25 Å 2 at 283 K) and (228.48 ± 0.27 Å 2 and 229.40 ± 0.23 Å 2 at 298 K), for the systems with DOD and SDS, respectively (Fig. 3 and Table S1). 3.2. Intermolecular interactions Due to the presence of many hydroxyl groups in the CDs and to their particular location, the most obvious specific contribution to the interaction between CD monomers and supramolecular structures with each other, as well as with water molecules, are the hydrogen bonds (H-bonds). Different contributions to the total number of H-bonds were calculated (see Methods section and Table S7). In general, the systems with bCD exhibit more H- bonds than those with a CD, clearly due to the larger number of hydroxyl groups of the former. Notheworthy, the number of H- bonds between CDs of different dimers of the same layer is marginal (Fig. 4), thus showing that this interaction is not the responsible for the formation and estabilization of the membranes, in contrast with what was suggested by other authors [4,8]. This happens for both CDs. Most hydrogen bonds are formed between the secondary hydroxyl groups of CDs of the same 2:1 complex -thus contributing to the stabilization of the supramolecular structures- and between the primary hydroxyl groups of the CDs and water molecules. Interestingly, the H-bonds between CDs of different monolayers are also marginal, the dominant interaction with the primary groups of the CDs being the H-bonds with water molecules, specially for the systems with SDS that systematically show more H-bonds with water than those with DOD. The total number of H-bonds between CDs and water molecules per complex is lower for the bilayer systems, regardless the ligand within the CDs. Altogether, the previous information suggests that an intermediate film of structured water could be mediating the interaction between 2:1 complexes of different monolayers. Density profile plots clearly show this water film (Fig. 4), which is stable throughout the whole trajectories involving 2:1 complex bilayers with both CDs. The behavior of the water molecules is highly sensitive to all the studied variables. For the simulations with a CD there is always a dry slice at the height-level corresponding to the interface between CD monomers forming the 2:1 complexes. This happens both in the systems with a single monolayer and in the systems with bilayers, which show one of these dry regions per layer (Figs. S7, S14, S22, S30, S37, S44, S52, S60). No waterrestricted regions are observed for any of the systems involving the wider bCD, although clear density-minima appear at the regions where both monomers are faced to each other in the complexes. This could be favored by the lower packing of the membranes based on this CD. The density of counterions for the simulations with SDS exhibits a maximum at the surface of the membranes, where the ionic heads are located. This happens on both sides of the membrane for the simulations with bilayers and just on the side corresponding to the SDS heads for the simulations with monolayers (Fig. 4 and Figs. S67, S74, S82, S90, S97, S104, S112, S120). The density of ions within the membrane is negligible for the systems with a CD and very low but not null for those with bCD. Fig. 3. Two views of the initial and final structures corresponding to SDS 1 CD 2 monolayers using a CD (left) and bCD (right) at 298 K. CDs that are closer and further from the SDS ionic head are represented in orange and yellow colors, respectively. The aliphatic chains of the SDS molecules are colored in blue, the sulfate group of SDS is in red, and sodium ions are shown as purple spheres. The central plot represents the average area per 2:1 complex (Inclusion Complex, IC) along the last microsecond of the sixteen MD trajectories studied. The different trajectories are indicated in the labels: monolayers (striped); bilayers (solid), systems at 283 K (dark color) and at 298 (light color). Fabián Suárez-Lestón, P.F. Garrido, Á. Piñeiro et al. Journal of Colloid and Interface Science 623 (2022) 938–946 941 3.3. Structure of the membranes Using the projection of the central point of each X 1 CD 2 complex on the XY plane as a reference (see methods section), Voronoi cells were determined for each frame of all the trajectories. The angle distribution for all the vertices was obtained and fitted to a series of skew normal functions to determine the percentage of hexagonal, square and rhombic lattice (see methods section, Fig. 5 and Table S6). Hexagonal lattices clearly dominate all the systems based on a CD, in agreement with the prediction of S. Yang et al. [8]. However, whereas these authors expected homogeneous rhombic geometries for the systems based on bCD, heterogeneous geometries were obtained in the present paper, probably due to the different conditions of the simulations compared to the experiments. The presence of square and rhombic cells is significant, in some cases comparable to the percentage of hexagonal lattice. In general, no clear pattern for the distribution between different types of geometries is observed for the systems based on bCD and the results do not depend on the guest molecule. This dynamic behaviour observed from the Voronoi cells of SDS 1 bCD 2 and DOD 1 - bCD 2 complexes contrasts with the extreme rigidity that have been associated with these systems by S. Yang et al. [8]. Our results suggest that the geometries induced by bCD in these structures look more like dynamic crystals instead of a crystal as a rigid ‘‘chemical cemetery” of molecules [33]. It can also be seen that for the systems with two layers of complexes based on a CD, the cyclodextrins of different layers are perfectly aligned (Figs. S16, S24, S46, S54). However, for the bilayers with bCD the complexes of each layer seem to be independent from those of the opposite one (Figs. S76, S84, S106, S114). Actually, no clear correlations were observed between structures of the two layers within the same trajectory. 3.4. Lateral and transversal diffusion Lateral diffusion coefficients of the complexes were obtained from the distribution of displacements in the XY plane over 1 ns and 5 ns time windows along the trajectories (see Fig. 6 and methods section). In general, the values obtained for the complexes in the monolayer systems were significantly lower (58–78% and 48–61% for time windows of 1 ns and 5 ns respectively) than those for the systems with bilayers (see Tables S2 and S3). This shows that, although the global structure of individual layers is independent from the presence of another one, there is a direct connection between them that is reflected in the movement of the molecules. The systems formed by bCD and DOD at 283 K and using a time window of 1 ns represent a singularity in this behavior since the difference between the diffusion coefficients of monolayers and bilayers is just 21%. In contrast, the same system at 298 K using a time window of 5 ns exhibits the opposite behavior, slowing down the displacement of complexes by 83% when a second layer is added. The replacement of DOD by SDS in the cavity of the CDs typically reduces the diffusion of the complexes. This change is typically moderate (4–10%) but in some cases is larger: bCD bilayers at 283 K for both time windows and bCD monolayers at 298 K for displacements along 5 ns. Notably, bCD monolayers at 283 K exhibit an increase in the diffusion coefficient upon the DOD-SDS replacement. The diffusion also depends clearly on temperature, as expected, typically increasing the coefficients by 40% when going from 283 K to 298 K. Again, the exception corresponds to the simulations with bCD and DOD which exhibit a significant larger (51%) and lower (15%) increase for monolayers using a time window of 1 ns and for bilayers using both time windows. The spontaneous diffusion of water molecules across the membranes based on a CD is much less favourable than in those with bCD. The total number of water molecules crossing the membrane during the 2- l s-long trajectories in the first case ranges from 200 water molecules for bilayers of DOD 1 a CD 2 complexes at 283 K to 1700 for monolayers of the same complex at 298 K. In contrast, between 20,000 and 300,000 water molecules completely cross the membranes based on bCD complexes (see Table S5). For the systems with a CD the diffusion of water molecules seems to be easier for monolayers with DOD 1 a CD 2 complexes than for those based on SDS 1 a CD 2 building blocks. However, this trend is reversed for bilayers of the same complexes. The behavior of the systems with bCD is similar to those with a CD but the number of water molecules crossing the membranes with the former CD is about 1000 times larger in all cases. This was already expected due to the presence of the water-forgiven region in the systems with a CD (see Figs. S7, S14, S22, S30, S37, S44, S52, S60). In this case, monolayers of complexes with both DOD and SDS allow passing more water molecules than bilayers of the same complexes and the water diffuses easier throughout systems with SDS than for those with DOD. As expected, the increase of temperature favours Fig. 4. (A) Different contributions to the mean number of H-bonds per CD: between monomers of the same 2:1 complex (red); between CDs of the same layer (yellow); and between CDs and water molecules (blue). The different trajectories are indicated in the labels: monolayers (striped); bilayers (solid), systems at 283 K (dark color) and at 298 (light color). (B) Water density profile as a function of time (greyscale) together with the average density profile of water molecules (orange), oxygen atoms joining the glucopyranoside groups (O4) of the CDs (green) and Na + ions (red); as a function of the simulation box Z coordinate, averaged over the last 1 l s of the trajectory corresponding to SDS 1 bCD 2 at 298 K. (C) Lateral view of the same membrane analyzed in (B) using the same representation code of Fig. 3. Fabián Suárez-Lestón, P.F. Garrido, Á. Piñeiro et al. Journal of Colloid and Interface Science 623 (2022) 938–946 942 the diffusion of water molecules throughout all membranes. Altogether, this information shows how these systems can be tuned to filter water, by changing temperature, the CD and the ligand. For the systems with SDS, the diffusion of counterions throughout the different membranes was also studied (Fig. 7, Figs. 132– 139 and Table S4). No crossing at all was observed for monolayers or bilayers based on a CD complexes, while the passing of ions in the systems with bCD complexes was significant. As in the case of water molecules, the temperature increase favours the diffusion of ions throughout the membranes. This happens both for monolayers and for bilayers. At 283 K, 47 events were observed for the monolayer system and only 22 for the bilayer while the difference in the number of ions crossing the monolayer and the bilayer was negligible at 298 K (71 and 70, respectively). This suggests how CD- based membranes could also be designed to filter specific ions. 4. Conclusions In agreement with previous experimental work [4,8], membranes based on DOD 1 bCD 2 and SDS 1 bCD 2 complexes in aqueous solution were stable for relatively long MD trajectories. Both, monolayers and bilayers of these supramolecular structures were tested here. Interestingly, similar membranes based on a CDs were also stable, with a more homogeneous lattice pattern but with completely different transport properties. Both monolayers and bilayers based on bCD are much more permeable to the passive diffusion of water and ions than those based on a CD. The interactions responsible for the stabilization of the membranes were also characterized. On the one hand, a specific pattern of H-bonds with many direct interactions between the secondary hydroxyls of CD-monomers of the same complex was observed. In contrast to Fig. 5. (A) Voronoi cells corresponding to the last frame of the simulation of SDS 1 a CD 2 (top) and SDS 1 bCD 2 (bottom) monolayers at 298 K. The shaded gray areas represent the projection of the unit cell simulation boxes on the XY plane. (B) Angles at the verticles of the Voronoi cells as a function of time together with the corresponding distributions for the same systems of (A). (C) Percentual contribution of hexagonal (yellow), square (purple), rhombic (green) and other types (white) of structures for each layer in the systems with a CD (left) and bCD (right). The mor bprefix of the system name indicates if the systems is a monolayer or a bilayer, respectively. The number between parenthesis is used to distinguish both layers in the bilayer systems. The colors indicate the temperature of the trajectory (blue for 283 K and red for 298 K). Fabián Suárez-Lestón, P.F. Garrido, Á. Piñeiro et al. Journal of Colloid and Interface Science 623 (2022) 938–946 943 what has been proposed as a driving force for the formation and stabilization of these membranes [4,8], just residual H-bonds were observed between CDs of different complexes. Given the presence of a pattern of highly ordered water molecules around CD complexes when they are isolated in aqueous solution [9], as well as the absence of H-bonds and also of water molecules between CDs of contiguous complexes in the studied membrane models, the most plausible interactions driving the formation and stabilization of these membranes is the favourable entropic contribution arising from the release of water molecules when the supramolecular complexes transfer from the bulk solution to the mesoscopic assemblies. The interaction with water molecules seems to be important for the stabilization of multilayers since the number of H-bonds between primary hydroxyl groups of different layers is residual and a layer of structured water molecules is present between both monolayers. This was observed for the bilayers involving both CDs and both ligands. X-ray scattering experiments of systems formed by SDS 1 bCD 2 complexes reveal a rhombic lattice structure with a = b = 1.52 nm (comparable to bCD diameter, 1.5 nm) and c = 104º[8]. The structure revealed by our analysis for membranes based on the same building block is a dynamic combination of hexagonal, square and rhombic lattices with a main cell parameter peak a = b = 1.6 1 nm, probably due to their fluid nature. Transitions between these different packing patterns were dynamically observed along the simulations (see Supplementary Movies), in agreemeng with the crystal adaptronics perspective [33]. For the systems with SDS 1 - a CD 2 complexes only the hexagonal pattern was observed for all the systems with almost identical lattice parameters a = b = 1.42, in good agreement with what had been predicted by S. Yang et al. [8]. Our results show that the passing of water and ions throughout the different membrane models can be tuned by changing the temperature, the CD and/or the ligand. CD replacement was found to be much more sensitive than the change of ligand or of temperature. The systems with a CD showed to be completely impermeable to ions, while they allow the crossing of a moderate amount of water molecules. In contrast, the systems with bCD are quite transparent to water and, to a lesser extent, to ions, although the diffusion can be modulated with temperature, with the addition of layers and with the replacement of the ligand forming the complexes. Only the passive diffusion of water and ions was considered in the present work but the differences are clear enough to get useful information. The natural continuation of this work requires filling an important experimental gap by synthesizing membranes based on a CD complexes. Simulations and experiments in the presence of external forces such as electric fields and different kind of ions would be useful to see the ability of these systems to filter specific charged particles. Overall, the information provided it the present work is expected to contribute to understand, at atomic resolution, the behavior of membrane systems formed by CD complexes, and thus Fig. 6. Lateral displacement distributions of DOD 1 a CD 2 inclusion complexes forming a monolayer at 283 K. The lines correspond to the fitting of the data calculated from the molecular dynamics trajectory over the last l s to the two-dimensional random walk equation (see methods section). Left and right plots correspond to calculations performed using 1 and 5 ns time windows, respectively. Fig. 7. Location of all Na + counterions (labelled by their index) in the simulation of the bilayer based on SDS 1 bCD 2 complexes at 298 K, as a function of time. White, grey and red regions correspond to the location of each ion above, below and within the membrane, respectively. Purple crosses indicate the initial time of an ion crossing. Direct transitions from grey to white, or vice versa, indicate the passing of ions throughout periodic boundaries of the simulation box. Fabián Suárez-Lestón, P.F. Garrido, Á. Piñeiro et al. Journal of Colloid and Interface Science 623 (2022) 938–946 944 to design new functional materials for specific applications based on them. CRediT authorship contribution statement F. S-L and P.F.G.: Formal analysis, Validation, Visualization. Á.P. and R G-F: Conceptualization, Funding acquisition, Project administration, Resources, Supervision. All authors: Investigation, Methodology, Writing - review & editing Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This work has received financial support from the Spanish Agencia Estatal de Investigación (AEI) and the European Regional Development Fund - ERDF (PID2019-111126RB-100, RTI2018- 098795-A-I00, and PID2019-111327GB-I00) and by the Xunta de Galicia (ED431F 2020/05 and Centro singular de investigación de Galicia accreditation 2019-2022, ED431G 2019/03) and the European Union (ERDF). R.G.-F. thanks Ministerio de Ciencia, Innovación y Universidades for a ‘‘Ramón y Cajal” contract (RYC- 2016-20335). P. F. G. thanks the Spanish Ministry of Economy and Competitiveness and the European Social Fund for his predoctoral research grant, reference BES-2016-076761. All calculations were carried out at the Centro de Supercomputación de Galicia (CESGA). Appendix A. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.jcis.2022.05.098. References [1] M.J. Baran, M.E. Carrington, S. Sahu, A. Baskin, J. Song, M.A. Baird, K.S. Han, K.T. Mueller, S.J. Teat, S.M. Meckler, C. Fu, D. Prendergast, B.A. Helms, Diversityoriented synthesis of polymer membranes with ion solvation cages, Nature 592 (7853) (2021) 225–231, https://doi.org/10.1038/s41586-021-03377-7. [2] W.J. Koros, C. Zhang, Materials for next-generation molecularly selective synthetic membranes, Nat. Mater. 16 (3) (2017) 289–297, https://doi.org/ 10.1038/nmat4805. [3] L. Jiang, Y. Peng, Y. Yan, M. Deng, Y. Wang, J. Huang, ‘‘Annular Ring” microtubes formed by SDS@2b-CD complexes in aqueous solution, Soft Matter 6 (2010) 1731–1736, https://doi.org/10.1039/B920608F. [4] L. Jiang, Y.u. Peng, Y. Yan, J. Huang, Aqueous self-assembly of SDS@2b-CD complexes: lamellae and vesicles, Soft Matter 7 (5) (2011) 1726–1731, https:// doi.org/10.1039/C0SM00917B. [5] B.G. Mathapa, V.N. Paunov, Cyclodextrin stabilised emulsions and cyclodextrinosomes, PCCP 15 (2013) 17903–17914, https://doi.org/10.1039/ C3CP52116H. [6] R.F. Fakhrullin, A.I. Zamaleeva, R.T. Minullina, S.A. Konnova, V.N. Paunov, Cyborg cells: Functionalisation of living cells with polymers and nanomaterials, Chem. Soc. Rev. 41 (2012) 4189–4206, https://doi.org/ 10.1039/C2CS15264A. [7] B.G. Mathapa, V.N. Paunov, Fabrication of viable cyborg cells with cyclodextrin functionality, Biomater. Sci. 2 (2) (2014) 212–219, https://doi.org/10.1039/ C3BM60162E. [8] S. Yang, Y. Yan, J. Huang, A.V. Petukhov, L.M.J. Kroon-Batenburg, M. Drechsler, C. Zhou, M. Tu, S. Granick, L. Jiang, Giant capsids from lattice self-assembly of cyclodextrin complexes, Nat. Commun. 2017 81 (8) (2017) 1–7, https://doi. org/10.1038/ncomms15856. [9] P.F. Garrido, M. Calvelo, R. Garcia-Fandiño, Á. Piñeiro, Rings, hexagons, petals, and dipolar moment sink-sources: The fanciful behavior of water around cyclodextrin complexes, Biomolecules 10 (2020) 431, https://doi.org/10.3390/ biom10030431. [10] P. Brocos, X. Banquy, N. Díaz-Vergara, S. Pérez-Casas, Á. Piñeiro, M. Costas, A critical approach to the thermodynamic characterization of inclusion complexes: Multiple-temperature isothermal titration calorimetric studies of native cyclodextrins with sodium dodecyl sulfate, J. Phys. Chem. B 115 (2011) 14381–14396, https://doi.org/10.1021/jp208740b. [11] P. Brocos, N. Díaz-Vergara, X. Banquy, S. Pérez-Casas, M. Costas, Á. Piñeiro, Similarities and differences between cyclodextrin-sodium dodecyl sulfate host-guest complexes of different stoichiometries: Molecular dynamics simulations at several temperatures, J. Phys. Chem. B 114 (2010) 12455– 12467, https://doi.org/10.1021/jp103223u. [12] H.J.C. Berendsen, J.P. Postma, W.F. van Gunsteren, J. Hermans, Interaction Models for Water in Relation to Protein HydrationIn, in: B. Pullman (Ed.), Intermol, Forces, Dordrecht, 1981, pp. 331–342, https://doi.org/10.1007/978- 94-015-7658-1_21. [13] N. Schmid, A.P. Eichenberger, A. Choutko, S. Riniker, M. Winger, A.E. Mark, W.F. van Gunsteren, Definition and testing of the GROMOS force-field versions 54A7 and 54B7, Eur. Biophys. J. 40 (7) (2011) 843–856, https://doi.org/ 10.1007/s00249-011-0700-9. [14] M.J. Abraham, T. Murtola, R. Schulz, S. Páll, J.C. Smith, B. Hess, E. Lindah, Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers, SoftwareX 1–2 (2015) 19–25, https://doi.org/10.1016/j.softx.2015.06.001. [15] H.J.C. Berendsen, D. van der Spoel, R. van Drunen, GROMACS: A messagepassing parallel molecular dynamics implementation, Comput. Phys. Commun. 91 (1-3) (1995) 43–56, https://doi.org/10.1016/0010-4655(95) 00042-E. [16] G. Bussi, D. Donadio, M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys. 126 (1) (2007) 014101, https://doi.org/10.1063/ 1.2408420. [17] M. Parrinello, A. Rahman, Polymorphic transitions in single crystals: A new molecular dynamics method, J. Appl. Phys. 52 (12) (1981) 7182–7190, https:// doi.org/10.1063/1.328693. [18] E. Mixcoha, J. Campos-Terán, Á. Piñeiro, Surface adsorption and bulk aggregation of cyclodextrins by computational molecular dynamics simulations as a function of temperature: a -CD vs b-CD, J. Phys. Chem. B 118 (25) (2014) 6999–7011, https://doi.org/10.1021/jp412533b. [19] Á. Piñeiro, J. Pipkin, V. Antle, R. Garcia-Fandino, Aggregation versus inclusion complexes to solubilize drugs with cyclodextrins. A case study using sulphobutylether-b-cyclodextrins and remdesivir, J. Mol. Liq. 343 (2021) 117588, https://doi.org/10.1016/j.molliq.2021.117588. [20] Á. Piñeiro, J. Pipkin, V. Antle, R. Garcia-Fandino, Remdesivir interactions with sulphobutylether-beta-cyclodextrin: a case study using selected substitution patterns, J. Mol. Liq. 346 (2022) 117157, https://doi.org/10.1016/ j.molliq.2021.117157. [21] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, Taylor & Francis Group, 1988, https://doi.org/10.1201/9780367806934. [22] T. Darden, D. York, L. Pedersen, Particle mesh Ewald: An Nlog(N) method for Ewald sums in large systems, J. Chem. Phys. 98 (12) (1993) 10089–10092, https://doi.org/10.1063/1.464397. [23] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, A smooth particle mesh Ewald method, J. Chem. Phys. 103 (19) (1995) 8577– 8593, https://doi.org/10.1063/1.470117. [24] S. Miyamoto, P.A. Kollman, SETTLE: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models, J. Comput. Chem. 13 (8) (1992) 952–962, https://doi.org/10.1002/jcc.540130805. [25] B. Hess, H. Bekker, H.J.C. Berendsen, J.G.E.M. Fraaije, LINCS: A linear constraint solver for molecular simulations, J. Comput. Chem. 18 (1997) 1463–1472, https://doi.org/10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0. CO;2-H. [26] R. Garcia-Fandiño, Á. Piñeiro, J.L. Trick, M.S.P. Sansom, Lipid Bilayer Membrane Perturbation by Embedded Nanopores: A Simulation Study, ACS Nano 10 (2016) 3693–3701, https://doi.org/10.1021/ACSNANO.6B00202/SUPPL_FILE/ NN6B00202_SI_001.PDF. [27] H. Leipold, E.A. Lazar, K.A. Brakke, D.J. Srolovitz, Statistical topology of perturbed two-dimensional lattices, J. Stat. Mech: Theory Exp. 2016 (4) (2016) 043103, https://doi.org/10.1088/1742-5468/2016/04/043103. [28] R.J. Gowers, M. Linke, J. Barnoud, T.J.E. Reddy, M.N. Melo, S.L. Seyler, J. Doman ´ski, D.L. Dotson, S. Buchoux, I.M. Kenney, O. Beckstein, MDAnalysis: A Python Package for the Rapid Analysis of Molecular Dynamics Simulations, in: Proc. 15th Python Sci. Conf., 2016, pp. 98–105, https://doi.org/10.25080/ MAJORA-629E541A-00E. [29] N. Michaud-Agrawal, E.J. Denning, T.B. Woolf, O. Beckstein, MDAnalysis: A toolkit for the analysis of molecular dynamics simulations, J. Comput. Chem. 32 (10) (2011) 2319–2327, https://doi.org/10.1002/jcc.21787. [30] C.R. Harris, K.J. Millman, S.J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N.J. Smith, R. Kern, M. Picus, S. Hoyer, M.H. van Kerkwijk, M. Brett, A. Haldane, J.F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, T.E. Oliphant, Array programming with NumPy, Nat. 585 (7825) (2020) 357– 362, https://doi.org/10.1038/s41586-020-2649-2. [31] J.D. Hunter, Matplotlib: A 2D graphics environment, Comput. Sci. Eng. 9 (3) (2007) 90–95, https://doi.org/10.1109/MCSE.2007.55. [32] P. Virtanen, R. Gommers, T.E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S.J. van der Walt, M. Brett, J. Wilson, K.J. Millman, N. Mayorov, A.R.J. Nelson, E. Jones, R. Kern, E. Larson, C.J. Carey, _ I. Polat, Y.u. Feng, E.W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E.A. Quintero, C.R. Harris, A.M. Archibald, A.H. Ribeiro, F. Pedregosa, P. van Mulbregt, A. Vijaykumar, A.P. Bardelli, A. Rothberg, A. Hilboll, A. Kloeckner, A. Scopatz, A. Lee, A. Rokem, C.N. Woods, C. Fulton, C. Fabián Suárez-Lestón, P.F. Garrido, Á. Piñeiro et al. Journal of Colloid and Interface Science 623 (2022) 938–946 945 Masson, C. Häggström, C. Fitzgerald, D.A. Nicholson, D.R. Hagen, D.V. Pasechnik, E. Olivetti, E. Martin, E. Wieser, F. Silva, F. Lenders, F. Wilhelm, G. Young, G.A. Price, G.-L. Ingold, G.E. Allen, G.R. Lee, H. Audren, I. Probst, J.P. Dietrich, J. Silterra, J.T. Webber, J. Slavic ˇ, J. Nothman, J. Buchner, J. Kulick, J.L. Schönberger, J.V. de Miranda Cardoso, J. Reimer, J. Harrington, J.L.C. Rodríguez, J. Nunez-Iglesias, J. Kuczynski, K. Tritz, M. Thoma, M. Newville, M. Kümmerer, M. Bolingbroke, M. Tartre, M. Pak, N.J. Smith, N. Nowaczyk, N. Shebanov, O. Pavlyk, P.A. Brodtkorb, P. Lee, R.T. McGibbon, R. Feldbauer, S. Lewis, S. Tygier, S. Sievert, S. Vigna, S. Peterson, S. More, T. Pudlik, T. Oshima, T.J. Pingel, T.P. Robitaille, T. Spura, T.R. Jones, T. Cera, T. Leslie, T. Zito, T. Krauss, U. Upadhyay, Y.O. Halchenko, Y. Vázquez-Baeza, SciPy 1.0: fundamental algorithms for scientific computing in Python, Nat. Methods 17 (3) (2020) 261–272. https:// www.nature.com/articles/s41592-019-0686-2. [33] P. Naumov, D.P. Karothu, E. Ahmed, L. Catalano, P. Commins, J. Mahmoud Halabi, M.B. Al-Handawi, L. Li, The Rise of the Dynamic Crystals, J. Am. Chem. Soc. 142 (31) (2020) 13256–13272, https://doi.org/10.1021/jacs.0c05440. Fabián Suárez-Lestón, P.F. Garrido, Á. Piñeiro et al. Journal of Colloid and Interface Science 623 (2022) 938–946 946