Cation–π interactions drive hydrophobic self-assembly and aggregation of niclosamide in water

Abstract

The beneficial medicinal effects of niclosamide have been reported to be hampered by poor aqueous solubility and so a higher concentration dosage is required. In this work, we have studied the aggregation properties of niclosamide in water by varying the number of monomers. We have employed all-atom classical molecular dynamics simulation in order to explore such properties. The equilibrium structure exists in an aggregated state with structural rearrangements of the stacking units. Niclosamide monomers tend to form clusters in an orderly manner and tend to aggregate in parallel and antiparallel orientations of the phenyl rings as the monomers are increased in number from 4 to 9. Upon increasing the size from 9 to 14, and from 49 to 150, a considerable dominance of the metastable parallel arrangement is observed, resulting in the formation of a closely packed cluster with hydrophobic contacts. The metastable conformation self-arranges to a T-shape before forming a stable planar antiparallel displaced conformation. The aggregated π–π parallel and cation–π antiparallel clusters in water exist in a β-conformer. We further observed that formation of a stable cluster aggregate entails the formation of an intermediate metastable cluster that disperses in solution forming a large stable cluster. We also discovered that movement of the water is faster in less aggregated clusters and as the cluster size increases, the mobility rate becomes much slower.

---

1 Introduction

Niclosamide (NC) (Fig. 1, 2,5-dichloro-4-nitrosalicylanilide) is a drug that the food and drug authority (FDA) has approved as an anthelminthic.^1,2 NC is included in the WHO’s list of essential medicines. It has so far been considered as one of the most highly potent multifunctional drugs.³ NC belongs to the Biopharmaceutical Classification System (BCS) Class II, meaning that it has a low solubility and high permeability.⁴ The crystalline form in which most drugs exist affects their solid-state characteristics.^5,6 Generally, the way drugs react plays a significant role on their performance characteristics. NC has been reported to crystallize into two solvate forms: two monohydrates (H_A and H_B) and one anhydrous.⁷ In the monohydrate form of NC, the water molecule is bound at the asymmetric unit of the drug surroundings by van der Waals envelopes of other drug molecules. The existence of NC in these different forms affects the dissolution characteristics and its bioavailability.^6,8

Fig. 1 (a) The chemical structure of niclosamide. (b) Crystal structure of two forms of niclosamide. The β-conformer is formed when the hydroxyl group (O20) from the ring forms a hydrogen bond with the hydrogen from a nitrogen (N11) at the middle, while the α conformer is formed when the hydroxyl group (O20) forms a hydrogen bond with the keto group (C12=O13). Atoms are shown by color: oxygen in red, nitrogen in blue, chlorine in green, carbon in grey and hydrogen in white.

In addition to its anthelminthic property, NC has been reported to be useful in the treatment of various categories of cancer and bacterial, fungal, and viral infections.^9–11 The efficacy of NC in the treatment of viral infections has been discussed in more detail, particularly in the treatment of coronavirus disease 19 (COVID-19) caused by SARS-CoV-2. The anti-viral properties of NC against SARS-CoV has been reported by Xu et al.¹² and Wu et al.,¹⁰ with an IC₅₀ of 1.56 μM. Also, NC was reported to inhibit the replication of MERS-CoV.¹³ Building from the existing knowledge about its potential anti-viral properties, NC was investigated for potential anti-viral activity against SARS-CoV-2.^14,15 Ko et al., reported the anti-SARS-CoV-2 activity of NC with an IC₅₀ of 0.28 μM.¹⁴ A recent study by Pindiprolu et al.¹⁶ has suggested possible mechanisms for the anti-SARS-CoV-2 activity of NC as blocking the endocytosis of SARS-CoV or preventing its autophagy through inhibition of the S-phase kinase-associated protein 2 (SKP2).¹⁶ The suggested mechanisms were developed from the analogy of MERS-CoV and SARS-CoV. Even though NC is a potent multifunctional drug, it is practically insoluble in water but slightly soluble in diethyl ether, ethanol, chloroform, tetrahydrofuran (THF) and dioxane.⁶ It has been reported that NC exhibits water solubilities of about 13.32 μg mL⁻¹ and 0.6–0.9 μg mL⁻¹ for the anhydrous and monohydrate forms, respectively.¹⁷ Due to its poor water solubility property the oral absorption of NC is limited. NC is mainly prepared in two types of dosage form: tablets and suspensions.⁷ The drug formulation of NC involves use of either NC anhydrate or the NC monohydrate forms.¹⁸ Due to the water affinity characteristic of the anhydrous crystal form of NC, the suspension develops into a cement-like sediment during formulation or when stored.¹⁷ Such limitations can be overcome by introducing an organic solvent, encapsulating with organic materials or coating with water-insoluble organic molecules. Various approaches have been applied to increase the solubility and hence the bioavailability of NC. They include: the use of lipid emulsions,¹⁹ solid–lipid nanoparticles,⁴ micelles,²⁰ nanoparticles,²¹ nanosuspensions,²² nanocrystals²³ and co-crystals.^18,24,25 However, little is known about how NC aggregates in an aqueous environment.

The aggregation process of NC needs to be addressed by atomistic simulation approaches which provide a better understanding of the characteristics of various materials at a molecular level. Computational investigation of biomolecular aggregation provides better microscopic details than the experimental approach. Problems related to aggregation and solubility have also been previously addressed using computational methods²⁶ for other systems. In order to address the problem of niclosamide aggregation at the atomistic level, we have carried out a molecular dynamics simulation of niclosamide monomers in water at an NPT ensemble for 50 ns.

2 Computational methods

2.1 Molecular dynamics simulation

In this study, molecular dynamics simulation was used to provide atomistic information on the hydrophobic interactions of niclosamide aggregation as well as the hydrophobic hydration that supports the stability. The structure of the niclosamide monomer was obtained from the pubChem database²⁷ in sdf format, energy minimized and converted to pdb file format. The coordinates and topology of niclosamide were obtained using the antechamber tool²⁸ together with the general AMBER force field (GAFF),²⁹ using the AM1-BCC charge method.³⁰ The obtained topology and coordinates were converted to Gromacs file format using the AMB2GMX python script. All the simulations reported were carried out at a temperature of 300 K and 1 bar pressure. Three individual systems comprising 4, 9 and 14 monomers of niclosamide were separately solvated in a cubic box containing a TIP4P water model.³¹ The box size was 3.0 × 3.1 × 3.0 nm for 4 monomers, 3.2 × 3.2 × 3.2 nm for 9 monomers, 3.5 × 3.5 × 3.5 nm for 14 monomers, 4.27 × 4.27 × 4.27 nm for 49 monomers and 10 × 10 × 10 nm for 150 monomers. The systems were energy minimized using a steepest descent method, and were equilibrated at an NPT ensemble for 500 ps. The temperature and pressure coupling were held using the v-rescale thermostat and Parinello–Rahman barostat, respectively.^32,33 A production run was carried out using the NPT ensemble for 50 ns with an integration time step of 2 fs. Periodic boundary conditions (PBCs) were applied in all directions. The Particle Mesh Ewald (PME) method was used to treat long-range electrostatic interactions with a cutoff distance of 11 Å for both electrostatic and van der Waals interactions, while covalent bonds were constrained using the LINCS algorithm.^34,35

The 2D free energy landscape was computed using the following equation;

F = −k_BT ln P_(r) (1)

where, k_B is the Boltzmann constant, T is the temperature and P_(r) is the probability distribution of the distances measured.

2.2 Quantum mechanical calculations

Geometry optimizations of the α and β-form were calculated using the Density Functional Theory (DFT) method at the M06-2X/def2-TZVP level of theory using the Gaussian 16 series of programs.³⁶ The Grimme’s dispersion correction³⁷ was also employed using the keyword empirical dispersion = gd3 in order to evaluate noncovalent interactions dominated by dispersion effects such as π-stacking. The interaction energies were corrected for the basis set superposition error (BSSE) using the counterpoise correction.³⁸ All the calculations were evaluated in an implicit water environment using the polarizable continuum model (PCM).³⁹ To account for solubility, the Gibbs free energy of solvation was computed at M062X/6-31G* using the implicit solvent (keyword “SCRF = SMD”).⁴⁰

3 Results and discussion

3.1 Hydrophobicity, niclosamide–niclosamide interaction and orientation

The radial distribution function (RDF) is an important tool for analysing the molecular interaction of niclosamide. Herein, the RDF was calculated by counting the number of atom pairs in niclosamide or niclosamide–water solvents at a specific distance. The density probability of atoms was then averaged as a function of distance. Fig. 2a shows the distribution of niclosamide molecules. It can be observed that the hydrophobic core of niclosamide with monomers arranged in parallel and antiparallel π–π stacking continuously reduces as the number of monomers increases from 4 to 14. The figure also shows three prominent peaks at three different positions suggesting the formation of three layers of clusters distributed in the system. Interestingly, we observed that, within the same system, the probability density peaks have different heights. This is attributed to the fact that as the cluster grows, it becomes more compact, and the monomers become closer to each other. However, as the system size increases from 9 to 14, the peak heights remain relatively similar suggesting that the system is closely packed (Fig. 2a).

Fig. 2 (a) Radial distribution function (RDF) for niclosamide systems (4 to 14), (b) dihedral angle distribution for niclosamide systems, the system formed by 14 niclosamide monomers is omitted for clarity. (c) Multi-step self-assembly process of niclosamide in water forming an aggregated state.

Upon comparing the three systems, one can observe that as the system size increases, the intensity of the RDF peak decreases. On the other hand, a high intensity RDF peak relates to the loose system with NC = 4, a similar observation has also been reported for a loose system of germanium atomic structures where the intensity of the RDF peak was also high.⁴¹

In order to understand how individual monomers adhere to each other in the aggregate, we calculated the dihedral angle distribution (C14, C12, N11 and C1, see Fig. 1a) that connects the two phenyl rings. Fig. 2b shows a bimodal distribution suggesting structural alignment of the two phenyl rings to each other. Visual inspection of the snapshot showed two types of phenyl–phenyl orientation existing in π–π antiparallel stacking and π–π parallel-displaced stacking (Fig. 2c). It is also important to highlight that, as the system increased from 4 to 9 the formation of parallel and antiparallel π–π stacking was equally observed from the equilibrium system (Fig. 2c). It is worth noting that, as the system size increased from 9 to 14, the metastable parallel conformation dominated, followed by the T-shaped conformation, which occurred for a short time before transforming to the stable anti-parallel conformation (Fig. S1†). The observed T-shape conformation clearly resembles those observed in proteins.⁴² The formation of this aggregate is facilitated by strong hydrophobic π–π interactions among the phenyl rings of niclosamide, providing increased stability.

Previous experimental studies on the crystal structure of niclosamide have reported two different results on conformational stability. Kosheleva et al.⁴³ crystalized niclosamide and observed two conformations, the α and β-conformation, and noted that the latter conformation was thermolabile and easily transformed to the stable α-conformation. Later on, Caira et al.⁸ crystalized niclosamide in different solvent and only observed the β-conformation. Our analysis based on the aggregated self-assembly behavior suggests the existence of the β-conformation only (Fig. 3a and b), which is in good agreement with experimental report by Caira et al.,⁸ who concluded that “the preferred conformation of the niclosamide molecule in the pseudopolymorphs investigated is the β-form”. It is worth noting that the pseudopolymorph packing of niclosamide in a different solvent reported by Caira et al.,⁸ is in good agreement with the packing aggregates reported in our work, and to some extent shows a similar packing motif to that observed in our simulation. The keto group formed hydrogen bonds with water molecules as seen in Fig. 3b, therefore stabilizing the observed π–π stacking, hence favouring the formation of the β-conformation only. These observations are similar to the experimental crystal structure of niclosamide in water, whereby the keto oxygen was observed to form a hydrogen bond with water.⁸ Although, in this study we have explored the self-assembly of niclosamide in water and extending and exploring the crystal structure stability at the atomistic level is warranted. This will further provide a better insight on the stability of the two forms of niclosamide in water.

Fig. 3 (a) A 2D free energy landscape showing the preferred conformation for niclosamide in water, (b) the intramolecular and intermolecular hydrogen bonds formed between the hydroxyl group and the nitrogen group of niclosamide as well as the keto group with water molecules. (c) Cluster size and transition formed during the aggregation process of niclosamide.

3.2 Cluster growth from segregation to stable aggregation

Cluster aggregation of several organic molecules such as caffeine has been studied using experimental methods. However, the transition process from segregation to aggregation is inadequately addressed from an experimental point of view. In this work, we have monitored the process of niclosamide aggregation using atomistic simulation.

We have also analyzed the cluster size and transition from one cluster to another to provide a better understanding. We observed a correlated dynamic behavior between cluster size and transitions; as the number of clusters decreased (from monomer to aggregated structure), the number of transitions also decreased, until a single cluster was formed and no further cluster growth was observed (Fig. 3c). The aggregation process of niclosamide resembles the nucleation process of a liquid, whereby each monomer tends to aggregate and segregate in the course of time until it attains an optimal cluster size.²⁶ A similar observation for curcumin resembling nucleation growth of liquid was previously studied using atomistic simulation.²⁶ The formation of stable aggregated niclosamide from the monomers was monitored by calculating the free energy landscape as a function of radius of gyration (R_g) and the root mean square deviation (RMSD) and solvent accessible surface area (SASA), which measures the degree of aggregation/segregation. Fig. 2c shows that the most aggregated system with 9 niclosamide monomers corresponds to lower R_g and RMSD values. However, this is a multi-step process involving the formation of three metastable steps before reaching the final cluster. The structures with larger R_g values are more segregated, however, with time, monomers assembled themselves until a final cluster was formed. The metastable states (Fig. 2c) showed that the monomers tended to join a larger stable state; while some monomers remained isolated for a longer time. Similarly, FES as the function of R_g and the SASA for the system with 14 niclosamide monomers shows that with large R_g values the system is in its segregated form. However, as the R_g values decreased the system attained larger stable aggregations. As observed from the 1D FES for the SASA (Fig. S2†), the most stable aggregate shows two minima with SASA values of 23 and 25 nm suggesting that the aggregated system might exist in two conformational states. Visual inspection of the aggregated form with SASA value = 25 nm (Fig. 4a) showed two conformational states, the parallel and anti-parallel, where the latter exists in a metastable state. The T-shaped conformation was also observed in the metastable state. Both parallel and T-shaped conformations were observed to be isolated and resulted in small aggregated clusters. We hypothesize that, after gaining some enthalpic and entropic stabilizations, the isolated monomers are able to bind to the large cluster and form a stable single aggregate. This is observed at a SASA value of 23 nm, where all monomers have joined, forming a stable anti-parallel conformation. In general, our MD calculations for 4, 9 and 14 monomers suggested that niclosamide forms a stable aggregated structure which tends to exist in an anti-parallel displaced form.

Fig. 4 (a) The formation of a stable aggregated system of niclosamide with 14 monomers. The aggregation is a multi-step process whereby individual monomers bind to form small aggregates dominated by the formation of parallel and T-shape conformations. With time, the small aggregates join to form a large cluster and re-orient to build a stable planar anti-parallel displaced conformation. (b) MEP surface of niclosamide, the energies shown at different points are in kcal mol⁻¹. (c and d) Interactions and arrangements of niclosamide dimers, the distance shown is in Å. To the right side of each figure is the NCI plot showing the interactions.

The systems discussed above helped to probe the interactions of the niclosamide molecule in water, however, in order to gain further insight into the aggregation behavior of niclosamide, large systems with 49 and 150 monomers of niclosamide were considered. Interestingly, as observed for systems with a small number of monomers, the larger system showed similar aggregation behavior and tendencies, however, with the formation of two well defined crystal like structures for the system with 49 monomers. Visual inspection showed that the aggregated systems are separated with bulk water, a phenomena which is not observed in small systems (Fig. 5a). To have a clear understanding of the interactions between niclosamide aggregates and water, density profiles were calculated along the z direction (Fig. 5b). Density profiles suggested that bulk water molecules tend to diffuse into the cluster surface of niclosamide over time. It is important to note from the density profile that, the two clusters formed are different sizes. As shown in the snapshot taken at 40 ns and in the density profile (Fig. 5a and b), it is evident that, the clusters formed do not have equal numbers of monomers, one cluster tends to have a larger number of well defined monomers with more penetrated hydrating water. The aggregation of 150 monomers further shows similar conformation (formation of antiparallel and parallel) arrangements with similar cluster growth observed compared to smaller systems, however, with a different nature. The larger system grows to form a chair like cluster, where small clusters join to form a larger thermodynamic cluster (Fig. 5).

Fig. 5 The aggregation behaviour of niclosamide in water. (a) The aggregation of 50 monomers of niclosamide, during the 60 ns niclosamide formed two clusters separated with water molecules. The two clusters are formed near the edge of the box leaving water molecules at the center. Bulk water molecules are also penetrating the large cluster to hydrate the other side. (b) Density profiles of niclosamide and water molecules at different times. (c) The aggregation process of 150 molecules of niclosamide, the process resembles the aggregation of 9 and 14 monomers, where small monomers joins to form a large cluster. During the 30 ns niclosamide monomers have formed a chair like cluster which continues to grow until 50 ns.

In order to confirm the observed orientational interactions, quantum mechanical calculations based on DFT were performed on various conformations. First, in order to know the electrostatic nature of the π-system and the most electron rich/deficient regions of β-niclosamide, the molecular electrostatic potential (MEP) was computed (Fig. 4b). As observed for the MD calculations, the plotted MEP surface shows a hydrogen bond formed between the NH–OH groups. Interestingly, we noted that the six membered ring connected to the nitroso group (NO₂) is less positive (+4.1 kcal mol⁻¹) than the other six membered ring (+5.7 kcal mol⁻¹) with the hydroxyl group (OH). This observation suggests that the formation of antiparallel cation–π or parallel π–π stacking displaced arrangements would be electrostatically favoured. As noted earlier, the antiparallel displaced dimer showed favourable interactions with the N atom of the nitroso group (which is less positive) which coincides with the location of the region of the more positive ring, forming a stable cation–π interaction. Eventually, this interaction is much favoured with a formation energy of −20.18 kcal mol⁻¹ compared to its counterpart, the non-displaced π–π stacking interaction, which resulted in an energy of −18.94 kcal mol⁻¹. However, the difference in energy is very minimal within the error bars. We further observed that the same conformation is likely to distort its planarity at the expense of forming an antiparallel hydrogen bond between the amide hydrogen and the carbonyl group of the other monomer. Such twisting and shifting reduces the previously observed cation–π interaction, which is the most stable, to the least stable π–π stacking interaction (Fig. 4c). A similar observation is noted for the parallel displaced conformation which is more stable (−15.40 kcal mol⁻¹) compared to the parallel conformation (−12.73 kcal mol⁻¹) (Fig. 4d). The observed antiparallel and parallel interactions were further characterized by computing the noncovalent interaction (NCI) plot index, which is an effective tool for showing noncovalent interactions. In Fig. 4c and d we show extended green isosurfaces located at the middle of the dimers characterized by acceptor and donor groups suggesting a large overlap of the π-interaction systems. It should be noted that the blue isosurfaces within the NCI plot represent formed intramolecular hydrogen bonds. Interestingly, our MD simulation is in reasonable agreement with the DFT calculation, where both suggest that the planar anti-parallel displaced conformation is the most stable. It should be noted that the observed T-shaped conformation is the metastable state changing from the parallel to anti-parallel conformation. Normally, the stable conformation for the niclosamide aggregate is stabilized by the formation of cation–π stacking interactions occurring in the displaced planar anti-parallel conformation. We have further observed that, on increasing the cluster size, the distances between planar fragments of parallel molecules decreased which indicates an indirect cooperativity in the hydrophobic contact (Fig. S4†). Experimental results on the crystal structure of the β-conformation showed that niclosamide exists in its planar conformation stabilized by the intramolecular hydrogen bond between the NH(11)–OH(20),⁸ which is in good agreement with both MD and DFT results.

3.3 Niclosamide cluster hydration

It has been noted that hydration has a great influence on cluster stability.^26,44 We examined structural water near aggregated niclosamide to gain a better understanding of the effect of hydration on cluster stability. Since niclosamide contains hydrophilic groups that tend to point outward, the cluster surface is hydrophilic. The observed local water structures in the vicinity of hydrophilic and hydrophobic surfaces exhibit different characteristics.

Fig. 6a, shows that solvated water molecules tend to penetrate into different levels of the cluster. Although the penetration of water molecules is limited to the hydrophobic region, they penetrate through the hydrophilic nitroso region to form a network of hydrogen bonds with the keto oxygen and the nitroso group as shown in Fig. 3b. The resulting hydrogen bond networks play an important role in stabilizing the cluster. To provide further understanding of the hydrogen bond lifetime, the hydrogen bond autocorrelation function, C(t) for continuous hydrogen bonding was calculated using eqn (2) and then fitted to the second order Legendre polynomial function as previously reported by Rawat and Biswas.⁴⁵ In eqn (2), S_i,j(t) measures the existence of hydrogen bonds between the ith and jth atoms at time t, and it will be zero if hydrogen bonds are absent. In Fig. 6b, we demonstrate that the intermolecular hydrogen bond relaxes faster for a small cluster i.e. with 9 niclosamide monomers, however; as the cluster size increases, the amount of penetrated water also increases thereby forming dense hydrogen bonds and hence slowing the water movement. The hydrophilic groups formed relatively more stable hydrogen bonds as compared to the hydrophobic region.

Fig. 6 (a) Hydration of the niclosamide cluster with interfacial water, (b) a hydrogen bond autocorrelation function for niclosamide aggregated clusters with 9 and 14 monomers of niclosamide. Hydrogen bonds persist for all the systems until 20 ps, then there are faster dynamics for a system with 9 monomers as compared to 14 monomers.

Fig. 8a, shows the water density near the cluster relative to the bulk solvent density. We observed that the cluster surfaces were hydrated. However, dewetting is observed to have taken place near the hydrophobic core. This is because hydrophobicity tends to break the already established hydrogen bond networks resulting in a decrease in density as the system size increases. It should also be noted that, although there is a wetting process, hydrophobic regions near the hydration site may affect the surrounding water density in the niclosamide cluster. We further observe that water in the hydration shell near the niclosamide cluster forms hydrogen bonds with the hydrophilic groups/sites. Fig. 8b shows the RDF between the niclosamide sites (C=O, O–H, –NH and –NO) and the water oxygen (OW). The only sites that form hydrogen bonds with the solvent are those belonging to the carboxylic and nitroso groups as seen in Fig. 3b. The carbonyl group shows a higher density of water indicating a strong hydrogen bond, followed by a hydroxyl group. The nitroso group (–NO) is found to have lower water density despite its exposure to the solvent. This is attributed to the reorientation of niclosamide groups. Furthermore, the –NH group was less hydrated and least exposed to the solvent surface as indicated in Fig. 8b.

3.4 Self-aggregation and solubility of niclosamide

The relationship between the self-association and the solubility issues of niclosamide in water is described. To gain insight into the solubility issues of aggregated niclosamide, we considered four different interaction modes observed in classical MD simulations of different cluster sizes, and calculated the Gibbs solvation free energy which relates to the solubility of niclosamide in water. The energy of interaction calculated from DFT for the four self-association modes; the antiparallel displaced (AP(dis)), antiparallel (AP), parallel (P) and parallel displaced (P(dis)) is shown in Fig. 4. As described above, the stable aggregate corresponds to the antiparallel displaced (−20.18 kcal mol⁻¹) and antiparallel (−18.94 kcal mol⁻¹) conformations. This suggests that the poor solubility of niclosamide in water is due to the antiparallel arrangements characterized by the formation of the cation–π and π–π interactions. To understand the relationship between the self-aggregation behavior and solubility issues, we have calculated the Gibbs free energy of solvation (ΔG_solv). Interestingly, the more stable aggregated form of niclosamide, that is, the antiparallel displaced mode (AP(dis)) exhibited the lowest (i.e. more positive value of) Gibbs free energy of solvation (−14.48 kcal mol⁻¹). On the other hand, the planar parallel aggregation mode showed increased solubility with a Gibbs solvation free energy of −22.47 kcal mol⁻¹ (Fig. 7). As mentioned earlier, niclosamide clusters are hydrated, Fig. 3b shows that the parallel cluster is the most hydrated and interacted with water compared to the antiparallel conformation, which is in good agreement with the calculated Gibbs free energy of solvation for the parallel clustered conformation. We suggest that, the marginal solubility (13.32 μg mL⁻¹) of niclosamide in water in its anhydrous form and 0.96 μg mL⁻¹ in its monohydrated forms is associated with the planar parallel structure and the formation of π–π interactions. On the other hand, the poor solubility is due to the antiparallel arrangement which results in the formation of cation–π interactions. We hypothesize that the aqueous solubility of niclosamide could be increased by breaking the cation–π interactions and favouring the formation of planar parallel interactions, however, this problem needs to be explored in future.

Fig. 7 The relationship between niclosamide self-aggregation and solubility in water. Stable aggregates exhibited the lowest Gibbs free energy of solvation and vice versa. AP(dis) means antiparallel displaced, AP means antiparallel, P means parallel and P(dis) means parallel displaced.

3.5 Water dynamics at niclosamide clusters

The hydrogen bond autocorrelation function provided an insight into the dynamical behavior of interfacial and bulk solvents. To further understand the dynamical behavior of interfacial water near the cluster interface, we calculated the residence time of water at the interface in the first hydration shell and the bulk. The residence time of water was estimated by calculating the survival probability N_w(t) as described by Rani et al.,⁴⁶ The residence time correlation function P_i(t_n, t) is the conditional probability for each ith water molecule; if P_i(t_n, t) = 1 it indicates that the water remains within the hydration shell in the time interval t_n and t_n + t, otherwise it is zero.

Fig. 8c and d shows the residence time correlation function of water at the cluster within the hydration shell and the bulk. The first hydration shell was defined for water within 3.5 Å and bulk water above 3.5 Å. Notably, different cluster sizes exhibited different dynamics of water. Although the hydrophobic tendency is pronounced, we observed that when niclosamide increased in size, the dynamics of water hydration reduced and remained intact for a longer time. At least 8 water molecules remained until 50 ps for the large cluster with 9 niclosamides. This observation is further supported by the calculated SASA values (Fig. S3†), where it is shown that the SASA values increase, while the cluster size increases. This implies that, an increase in size resulted in an increase of trapped water molecules which are responsible for hydrating the surface. This observation is similar to the hydrogen bond autocorrelation function. However, the dynamical behaviour of the bulk solvent did not show much difference for different clusters (Fig. 8c).

Fig. 8 (a) The radial distribution function (RDF) between niclosamide and water oxygen (OW) for different systems, (b) the RDF between niclosamide sites (C=O, –OH, –NH, and –NO) and water oxygen (OW), (c) dynamics of the hydration water molecules that remained trapped within the first shell (3.5 Å), (d) dynamics of bulk water.

To gain insight into the timescale associated with the dynamics of water, the residence time correlation function was fitted to an exponential function as described previously by Rani et al.,⁴⁶

S(t) = a exp(−(t/τ_s)^γ) + b exp(−(t/τ₂)) + c exp(−(t/τ₃)) + n_p (4)

where, n_p is the number of water molecules present at niclosamide throughout the simulation time, τ_s is the residence time for the decay of the stretched exponential. τ₂ and τ₃ are the residence times for slow biexponential decay for the first and second components, respectively. γ represents a quantitative measure of the deviation of the relaxation curve from a classical exponential function. The value of γ = 1 for a classical exponential function. A large γ value indicates the presence of system temporal disorders.⁴⁶ The best fit parameters are presented in Table 1, the average residence time of the stretched exponential is obtained using the expression:⁴⁶where, Γ is the gamma function.

Table 1 Residence time of interfacial and bulk water calculated by fitting the survival probability

Cluster	a	〈τs〉	γ	b	τ2	c	τ3	np
NC no = 4 (3.5 Å)	14.2502	17.3738	1.42502	5.36237	3.98729	52.1834	8.89985	0.00068
NC no = 9 (3.5 Å)	23.2337	20.3232	2.32337	17.5998	6.15648	58.7136	10.3767	0.00044
NC no = 14 (3.5 Å)	19.5257	27.0507	1.95273	20.9493	5.94705	99.1797	12.1214	0.00333
NC no = 9 (bulk)	40.0821	25.9479	4.00821	6.59938	−0.417845	65.3204	12.0326	0.01359

---

From Table 1, it can be observed that different cluster sizes exhibited different γ values. This is expected in light of the differences in the cluster sizes and molecule behavior. Table 1 shows that the system with a large number of niclosamide monomers exhibited the slowest dynamics of water, whilst for a system with few niclosamide monomers it exhibited the fastest dynamics. The general trend is that, as the cluster size increases, the ability to trap water molecules inside the shell increases. The bulk water exhibits a different dynamical behavior as compared to the interfacial water molecules near the cluster.

4 Conclusions

How water mediates the interaction/assembly of a biological system is an intriguing question that both experimentalists and computational theorists have yet to address precisely. To address this question at the molecular level, we have simulated niclosamide aggregation in water, maybe for the first time. Niclosamide tends to form a large cluster in an orderly manner. Several intermediate metastable clusters are formed during the development of stable large clusters. The dynamical transitions from segregation to aggregation are also monitored. We have also studied the dynamics of water hydration at the interface of the clusters. We showed that water molecules at the interface present rather different dynamical properties compared to the bulk water. The clusters formed have a defined crystalline form. Although this work has described the hydrophobic aggregation of niclosamide in water, problems relating to heterogeneous crystal water interaction, polymorphism growth, solubility measurements, and how they can be addressed by quantitative approaches like the potential mean of force for large hydrophobic clusters are worth pursuing in future.

Author contributions

All authors conceived and designed the study, DMS and LWK performed MD and DFT experiments and DMS drafted the manuscript. All authors read, revised and approved for submission.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors would like to thank the South Africa Centre for High Performance Computing (CHPC) for providing computational resources on their Lengau cluster for this research. This study was supported by the University of Dodoma through the SAS funding program.

References

1. P. Andrews, J. Thyssen and D. Lorke, Pharmacol. Ther., 1982, 19, 245–295.
2. R. D. Pearson and E. L. Hewlett, Niclosamide therapy for tapeworm infections, 1985.
3. W. Chen, R. A. Mook Jr, R. T. Premont and J. Wang, Cell. Signalling, 2018, 41, 89–96.
4. M. U. Rehman, M. A. Khan, W. S. Khan, M. Shafique and M. Khan, Artif. Cells, Nanomed., Biotechnol., 2018, 46, 1926–1934.
5. A. Burger and R. Ramberger, Microchim. Acta, 1979, 72, 259–271.
6. R. V. Manek and W. M. Kolling, AAPS PharmSciTech, 2004, 5, 101–108.
7. M. M. de Villiers, M. D. Mahlatji, E. C. van Tonder, S. F. Malan, A. P. Lötter and W. Liebenberg, Drug Dev. Ind. Pharm., 2004, 30, 581–592.
8. M. Caira, E. Van Tonder, M. De Villiers and A. Lötter, J. Inclusion Phenom. Mol. Recognit. Chem., 1998, 31, 1–16.
9. Z. Lou, Y. Sun and Z. Rao, Trends Pharmacol. Sci., 2014, 35, 86–102.
10. C.-J. Wu, J.-T. Jan, C.-M. Chen, H.-P. Hsieh, D.-R. Hwang, H.-W. Liu, C.-Y. Liu, H.-W. Huang, S.-C. Chen and C.-F. Hong, et al., Antimicrob. Agents Chemother., 2004, 48, 2693–2696.
11. A. Jurgeit, R. McDowell, S. Moese, E. Meldrum, R. Schwendener and U. F. Greber, PLoS Pathog., 2012, 8, e1002976.
12. J. Xu, P.-Y. Shi, H. Li and J. Zhou, ACS Infect. Dis., 2020, 6, 909–915.
13. N. C. Gassen, D. Niemeyer, D. Muth, V. M. Corman, S. Martinelli, A. Gassen, K. Hafner, J. Papies, K. Mösbauer and A. Zellner, et al., Nat. Commun., 2019, 10, 1–16.
14. M. Ko, S. Y. Chang, S. Y. Byun, A. Ianevski, I. Choi, P. H. d. d’Orengiani, E. Ravlo, W. Wang, M. Bjørås and D. E. Kainov, et al., Viruses, 2021, 13, 651.
15. G. Wang, H. Gaikwad, M. McCarthy, M. Gonzalez-Juarrero, Y. Li, M. Armstrong, N. Reisdorph, T. Morrison and D. Simberg, Precis. Nanomed., 2021, 4, 724–737.
16. S. K. S. Pindiprolu and S. H. Pindiprolu, Med. Hypotheses, 2020, 140, 109765.
17. E. C. van Tonder, T. S. Maleka, W. Liebenberg, M. Song, D. E. Wurster and M. M. de Villiers, Int. J. Pharm., 2004, 269, 417–432.
18. P. Sanphui, S. S. Kumar and A. Nangia, Cryst. Growth Des., 2012, 12, 4588–4599.
19. X. Zhang, Y. Zhang, T. Zhang, J. Zhang and B. Wu, J. Microencapsulation, 2015, 32, 496–502.
20. A. Russo, D. S. Pellosi, V. Pagliara, M. R. Milone, B. Pucci, W. Caetano, N. Hioka, A. Budillon, F. Ungaro and G. Russo, et al., Int. J. Pharm., 2016, 511, 127–139.
21. S. Naqvi, S. Mohiyuddin and P. Gopinath, R. Soc. Open Sci., 2017, 4, 170611.
22. G. Costabile, I. d’Angelo, G. Rampioni, R. Bondì, B. Pompili, F. Ascenzioni, E. Mitidieri, R. d’Emmanuele di Villa Bianca, R. Sorrentino and A. Miro, et al., Mol. Pharmaceutics, 2015, 12, 2604–2617.
23. Y. Ye, X. Zhang, T. Zhang, H. Wang and B. Wu, Drug Dev. Ind. Pharm., 2015, 41, 1416–1424.
24. F. Grifasi, M. R. Chierotti, K. Gaglioti, R. Gobetto, L. Maini, D. Braga, E. Dichiarante and M. Curzi, Cryst. Growth Des., 2015, 15, 1939–1948.
25. D. Luedeker, R. Gossmann, K. Langer and G. Brunklaus, Cryst. Growth Des., 2016, 16, 3087–3100.
26. M. K. Hazra, S. Roy and B. Bagchi, J. Chem. Phys., 2014, 141, 18C501.
27. S. Kim, J. Chen, T. Cheng, A. Gindulyte, J. He, S. He, Q. Li, B. A. Shoemaker, P. A. Thiessen and B. Yu, et al., Nucleic Acids Res., 2019, 47, D1102–D1109.
28. A. W. S. Da Silva and W. F. Vranken, BMC Res. Notes, 2012, 5, 1–8.
29. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
30. A. Jakalian, D. B. Jack and C. I. Bayly, J. Comput. Chem., 2002, 23, 1623–1641.
31. M. A. González and J. L. Abascal, J. Chem. Phys., 2011, 135, 224516.
32. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
33. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
34. T. I. Cheatham, J. Miller, T. Fox, T. Darden and P. Kollman, J. Am. Chem. Soc., 1995, 117, 4193–4194.
35. B. Hess, H. Bekker, H. J. Berendsen and J. G. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
36. M. Frisch, G. Trucks, H. Schlegel, G. H. Scuseria, M. Robb, J. Cheeseman, G. Scalmani, V. Barone, G. H. Petersson, Nakatsuji, et al., Gaussian, Inc., Wallingford CT, 2016.
37. A.-R. Allouche, J. Comput. Chem., 2011, 32, 174–182.
38. S. F. Boys and F. Bernardi, Mol. Phys., 1970, 19, 553–566.
39. E. Cances, B. Mennucci and J. Tomasi, J. Chem. Phys., 1997, 107, 3032–3041.
40. A. V. Marenich, C. J. Cramer and D. G. Truhlar, J. Phys. Chem. B, 2009, 113, 6378–6396.
41. P. Sneha and C. G. P. Doss, Adv. Protein Chem. Struct. Biol., 2016, 102, 181–224.
42. G. B. McGaughey, M. Gagné and A. K. Rappé, J. Biol. Chem., 1998, 273, 15458–15463.
43. L. Kosheleva and A. Bekhli, Pharm. Chem. J., 1974, 8, 576–578.
44. M. Vazdar, J. Heyda, P. E. Mason, G. Tesei, C. Allolio, M. Lund and P. Jungwirth, Acc. Chem. Res., 2018, 51, 1455–1464.
45. N. Rawat and P. Biswas, J. Phys. Chem. B, 2014, 118, 3018–3025.
46. P. Rani and P. Biswas, J. Phys. Chem. B, 2015, 119, 10858–10867.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1ra05358b

---