Effect of Choice of Solvent on Crystallization Pathway of Paracetamol: An Experimental and Theoretical Case Study

The choice of solvents influences crystalline solid formed during the crystallization of active pharmaceutical ingredients (API). The underlying effects are not always well understood because of the complexity of the systems. Theoretical models are often insufficient to describe this phenomenon. In this study, the crystallization behavior of the model drug paracetamol in different solvents was studied based on experimental and molecular dynamics data. The crystallization process was followed in situ using time-resolved Raman spectroscopy. Molecular dynamics with simulated annealing algorithm was used for an atomistic understanding of the underlying processes. The experimental and theoretical data indicate that paracetamol molecules adopt a particular geometry in a given solvent predefining the crystallization of certain polymorphs.


Introduction
Polymorphs are crystalline modifications exhibiting the same composition, but different packing arrangements. Polymorphism is a universal phenomenon that can be found in metals, minerals, and organic molecules. A prominent example of organic polymorphism is cocoa butter, a key ingredient in chocolate. At least six different crystal forms are known differing in properties like appearance, texture, and melting point [1,2]. Polymorphism is also the reason behind a phenomenon where a desirable form in producing a drug could be no longer produced in the same facility, known as disappearing polymorphs [3,4]. One well-known example in this context is Ritonavir, a drug used in the treatment of HIV infection, failed the dissolution specifications test because of the presence of a polymorph having approximately 50% lower intrinsic solubility than the reference form [5]. This example shows that crystallization processes and the entanglement of physicochemical interactions need to be understood to decipher the fundamental formation processes. Crystallization yields one of the crystal forms out of many polymorphs depending on various parameters such as cooling rate or choice of solvent [6]. A crystal may be considered as the embodiment of a supra-molecular assembly whose synthons are the molecules [7]. Davey et al. found that the self-associates in solution bear a resemblance to the structural synthons of crystal structures in a number of organic systems [8]. The crystallization process involves two steps: nucleation and crystal growth, out of which nucleation is the rate-determining step [9]. Different theories postulating different pathways have been proposed explaining the crystallization process [10]. Classical nucleation theory (CNT) and multinucleation pathways are discussed in this context [11][12][13][14][15][16][17]. It has been proven that a single system of solute-solvent can behave differently and The crystallographic structure of paracetamol shows an extensive hydrogen-bonding network, contributing to the stability of the structure [33][34][35][36][37]. Due to this polymorphism, the monoclinic form I and orthorhombic form II are processed differently during the formulation of the paracetamol tablets. The form I requires extra binder during punching of the tablets, because of its corrugated packing, usually resulting in poor compression properties. The form II can be directly compressed into a tablet, which has a parallel hydrogen-bonded sheets-like structure [38]. It was reported that the transition of paracetamol from form II to form I did not affect bioavailability [39].
In this work, we investigate the effects of solvent selection on the crystallization of paracetamol from the solution, using solvent mixtures and dilutions. An ultrasonic trap was used as sample environment to ensure minimal extrinsic influences. The significant similarity in structural synthon correspondence was further revealed by in situ solution spectroscopy and computer simulation studies. The experimental data is backed by molecular dynamics data. We found that the choice of solvent is a driving force towards a particular arrangement of paracetamol molecules triggering the crystallization of a certain polymorph. Furthermore, the solvents exhibit different strengths to drive crystallization towards a certain endpoint in case of paracetamol.

Crystallization Experiments
Crystallization studies from liquid solutions were performed in a custom-made acoustic levitator. The acoustic levitator allows conducting contact-free crystallization studies and in situ measurements. A droplet of the solution can be fixed in a stable and undisturbed position by means of an ultrasonic field. The environment around the sample can be controlled regarding the surface, temperature, and humidity by passing a cold/hot stream of nitrogen [40]. This setup allows studying the phenomena of crystallization, minimizing influencing factors such as extrinsic heterogeneous seeding. During the experiment, the solvent evaporated and led to a gradual increase in the concentration of the droplet, which finally crystallized. The custom-made acoustic levitator is shown schematically in Figure 1. The droplet and eventually the crystallized paracetamol was held in the trap using a sonotrode vibrating at a frequency of 58 kHz and a concave reflector, which created a standing acoustic wavefield with a defined number of local nodes. Droplets or solid crystals with diameters less than half the wavelength can be levitated just below in either of the nodes of the wave due to axial radiation pressure and radial Bernoulli stress.
Crystals 2020, 10, x FOR PEER REVIEW 3 of 11 environment to ensure minimal extrinsic influences. The significant similarity in structural synthon correspondence was further revealed by in situ solution spectroscopy and computer simulation studies. The experimental data is backed by molecular dynamics data. We found that the choice of solvent is a driving force towards a particular arrangement of paracetamol molecules triggering the crystallization of a certain polymorph. Furthermore, the solvents exhibit different strengths to drive crystallization towards a certain endpoint in case of paracetamol.

Crystallization Experiments
Crystallization studies from liquid solutions were performed in a custom-made acoustic levitator. The acoustic levitator allows conducting contact-free crystallization studies and in situ measurements. A droplet of the solution can be fixed in a stable and undisturbed position by means of an ultrasonic field. The environment around the sample can be controlled regarding the surface, temperature, and humidity by passing a cold/hot stream of nitrogen [40]. This setup allows studying the phenomena of crystallization, minimizing influencing factors such as extrinsic heterogeneous seeding. During the experiment, the solvent evaporated and led to a gradual increase in the concentration of the droplet, which finally crystallized. The custom-made acoustic levitator is shown schematically in Figure 1. The droplet and eventually the crystallized paracetamol was held in the trap using a sonotrode vibrating at a frequency of 58 kHz and a concave reflector, which created a standing acoustic wavefield with a defined number of local nodes. Droplets or solid crystals with diameters less than half the wavelength can be levitated just below in either of the nodes of the wave due to axial radiation pressure and radial Bernoulli stress. The climate unit provides constant relative humidity and temperature. In conjunction with the acoustic levitator, time-resolved Raman spectroscopy was used. The Raman measurements were done using Raman RXN1 TM instrument (Kaiser optical systems) with an excitation wavelength of 785 nm. The humidity RH and temperature T were maintained at 17.5 ± 2.5% and 22.0 ±1.0 °C during the measurements.
Materials: paracetamol (>99% purity, Sigma-Aldrich, St. Louis, US) was used without further purification. X-ray diffraction and Raman spectroscopy were used for the characterization of the initial crystalline structure, which was found out to be the monoclinic form I. Methanol (>99.99%) was procured from Carl Roth GmbH + Co. KG (Karlsruhe, Germany), ethanol (>99.8%) from Th. Geyer GmbH + Co. KG (Renningen, Germany), and chloroform (>99.99%) from AppliChem GmbH (Darmstadt, Germany) The climate unit provides constant relative humidity and temperature. In conjunction with the acoustic levitator, time-resolved Raman spectroscopy was used. The Raman measurements were done using Raman RXN1 TM instrument (Kaiser optical systems) with an excitation wavelength of 785 nm. The humidity RH and temperature T were maintained at 17.5 ± 2.5% and 22.0 ± 1.0 • C during the measurements.

Conformer Preparation
The Cambridge Structural Database (CSD) entries HXACAN09 and HXACAN23 were chosen as representatives for the monoclinic (form I) and orthorhombic (form II) paracetamol polymorphs, respectively [41,42]. Single-molecule structures were prepared. For multiple copies of the paracetamol structures during the simulations, genbox algorithm from GROMACS was used, which allows us to randomly add the structures in the simulation box of a defined size. Polymorph I and II of paracetamol have a unique arrangement of molecules within their crystal structure. We chose a synthon within each of these structures, which are explicit structural identifiers of polymorphic crystal structures, and exclusive in their relative crystal structures. In form I, (CSD entry: HXACAN09, Figure 2a) paracetamol molecules forming hydrogen bonds were arranged orthogonal to each other.

Conformer Preparation
The Cambridge Structural Database (CSD) entries HXACAN09 and HXACAN23 were chosen as representatives for the monoclinic (form I) and orthorhombic (form II) paracetamol polymorphs, respectively [41,42]. Single-molecule structures were prepared. For multiple copies of the paracetamol structures during the simulations, genbox algorithm from GROMACS was used, which allows us to randomly add the structures in the simulation box of a defined size. Polymorph I and II of paracetamol have a unique arrangement of molecules within their crystal structure. We chose a synthon within each of these structures, which are explicit structural identifiers of polymorphic crystal structures, and exclusive in their relative crystal structures. In form I, (CSD entry: HXACAN09, Figure 2a) paracetamol molecules forming hydrogen bonds were arranged orthogonal to each other. A vector (black arrow) passing from the nitrogen atom at one terminal of the aromatic ring towards oxygen atom of the hydroxyl group at the opposite side of the ring cut the similar vector from another molecule in an orthogonal fashion. Contrarily, form II (CSD entry: HXACAN23, Figure  2, b) exhibited hydrogen bonding in a completely different manner. The vectors drawn in similar fashion from nitrogen atom towards oxygen atom attached to ring from of two molecules forming hydrogen bonds in the crystal structure were parallel to each other. This fundamental difference in the crystal structure of two polymorphs was exploited in the presented study.

Simulation Protocol
The aggregation of molecules via non-covalent interaction is a process of fundamental importance for crystallization. During such a process, energy is released, a local or global optimum is achieved by the system, which dictates the nature of the polymorph. In mathematical equivalence, such a process could be rendered as an optimization problem. There are a number of ways this problem can be tried to be solved, and simulated annealing is one of the widely used algorithms used for such an optimization purpose. It can be excellently implemented for molecular systems. This process is akin to annealing in metallurgy, where metal is heated and cooled in a controlled manner to increase the size of the crystals and reduce defects. Enough time is spent at the cooling cycle to reach thermal equilibrium in order for molecules to reach a thermal equilibrium where energy levels A vector (black arrow) passing from the nitrogen atom at one terminal of the aromatic ring towards oxygen atom of the hydroxyl group at the opposite side of the ring cut the similar vector from another molecule in an orthogonal fashion. Contrarily, form II (CSD entry: HXACAN23, Figure 2b) exhibited hydrogen bonding in a completely different manner. The vectors drawn in similar fashion from nitrogen atom towards oxygen atom attached to ring from of two molecules forming hydrogen bonds in the crystal structure were parallel to each other. This fundamental difference in the crystal structure of two polymorphs was exploited in the presented study.

Simulation Protocol
The aggregation of molecules via non-covalent interaction is a process of fundamental importance for crystallization. During such a process, energy is released, a local or global optimum is achieved by the system, which dictates the nature of the polymorph. In mathematical equivalence, such a process could be rendered as an optimization problem. There are a number of ways this problem can be tried to be solved, and simulated annealing is one of the widely used algorithms used for such an optimization purpose. It can be excellently implemented for molecular systems. This process is akin to annealing in metallurgy, where metal is heated and cooled in a controlled manner to increase the size of the crystals and reduce defects. Enough time is spent at the cooling cycle to reach thermal equilibrium in order for molecules to reach a thermal equilibrium where energy levels follow Boltzmann distribution. This view was adopted by some authors, but it was strongly argued against in molecular dynamics situations, where equilibrium may not be achieved or even desired in a system like ours. We did not apply simulated annealing in the classical sense. This means, we did not expect an equilibration of the system and tried to extract canonical ensemble statistics from it. Instead, we conjectured, that due to the "geometry" of the solvent, paracetamol molecules can only "meet" in specific arrangements (angles) that vary between the solvents. Thus, instead of waiting "till equilibration", we waited "till synthon forming" and then listed the different possible angles of the pairs of paracetamol molecules.
The goal here is to gather adequate statistical weightage. We used GROMACS 4.6 for simulating paracetamol molecules in different solvents. The paracetamol and solvent molecules were parameterized using amber99sb forcefield, which is optimized and suitable for organic molecules containing C, H, O, and N elements. Sixteen paracetamol molecules were added along with 1528 ethanol molecules, 1379 methanol molecules, 871 chloroform molecules, and 5337 water molecules. Although this addition is random, general numbers decrease with a relative increase in the molecular size of the solvent. Exact molar concentrations of solutes and solvents in the experimental setup were sought to match in simulations; however, such attempts were not successful because of overfilling in the simulation box. Octahedron shaped simulation box of the longest internal diagonal of 6 nanometers with periodic boundary conditions was set up for simulations. For the energy minimization, the conjugate gradient method was used with 15,000 steps without any constrains.
Once energy minimization ensured that we had a reasonable starting structure, for equilibration of the system, leap-frog integrator was used with each time step lasting for 100 ps at 293 K. NPT ensemble, also known as Isothermal-Isobaric ensemble, where a number of particles applied pressure and temperature at the given time remain constant, was used for this purpose. The trajectory was generated using simulated annealing cycles starting from 100 to 600 K (refer Table 2) for 110 ps, and were defined for solute and solvent. A total of 45 cycles of simulated annealing were simulated. This procedure was intended for rapid, periodic heating and cooling of the system. Simulated annealing helps to sample as many as possible conformations of paracetamol and solvent molecules, thus allowing more conformational coverage. As a result, these types of simulations offered a good statistical sampling. Subjecting the system to periodic high temperatures also destabilizes intermolecular hydrogen bonds, which are formed in cycles before, hence the stagnancy in the system was forbidden, and new combinations of interactions during the upcoming cycles in the simulation were made possible. The system was simulated for 5 nanoseconds (simulation time). The trajectory was processed to get xtc output file, with nojump argument in order to keep molecules intact. This helps to avoid fragmentation of molecules at periodic boundaries, to make sure there are no errors in the analysis because of the formation of new species. A separate MATLAB parser code named as Gro2mat was used to convert the .xtc file to easy xyz co-ordinates, which were readable in MATLAB [43]. Co-ordinates of paracetamol molecules were skimmed off from each step of the output file and were further processed.

Analysis of the Trajectory
After the MD simulation, the generated trajectories of molecular states were analyzed. The effect of the solvent on the arrangement of the solute molecules was the desired quality of interest. To quantify the strength of this effect, we investigated whether the paracetamol molecules were organized according to the synthons identified in form I or to form II. This question can only be answered if two paracetamol molecules approach each other and form a hydrogen bond because the stability of the form I and form II is due to a special pattern of hydrogen bonds. For this analysis, we defined a set Crystals 2020, 10, 1107 6 of 10 of possible hydrogen-bond donor and hydrogen-bond acceptor atoms of the paracetamol molecules within the simulation. Whenever a donor and an acceptor atom approach each other less than 2.5 Å during simulation, we probed whether the current arrangement of these two respective molecules corresponded to form I or form II. Nitrogen atom bound to the C6-ring and the opposite oxygen atom bound to the C6-ring defined an orientation vector in each molecule. Whenever an H-bond between two molecules was logged, we measured the angle between the two orientation vectors spanned by the described atoms. If this angle was between 45 and 135 • , which is half the range of the whole analytical angle. We then denoted the corresponding arrangement as synthon of form I. Otherwise, we denoted it as synthon of form II. For details of the method, refer the Supplementary Material.
The simulation studies carried out of paracetamol molecules in the solvents previously showed that the varying strength (in terms of numbers) of the hydrogen bonds could indicate the intrinsic ability of solvents to drive the crystallization of paracetamol towards a particular polymorph. To explore this phenomenon in detail, an empirical serial dilution series was created, and paracetamol could crystallize out of these mixed solvents.

Results
Our experimental findings concurred with findings of previous studies of paracetamol polymorphs in specific solvents listed in Table 3. Table 3. Polymorphs of paracetamol obtained from different solvents (experimental).

Ethanol Methanol Chloroform Water
Results of MD simulations of paracetamol in pure solvents indicated two elemental effects: the arrangement of paracetamol molecules with respect to each other and second bound by hydrogen bond; and the number of hydrogen bonds in each of the solvents. Every time a hydrogen bond was detected in the simulation, the stated angle between vectors of participating paracetamol molecules was recorded. This process was indicative of searching synthons in the trajectory produced in MD simulations. The number of hydrogen bonds formed within these synthons in each of these solvents varied significantly. These are tabulated in Table 4. The angles between a pair of paracetamol molecules were then represented in the form of histograms for different solvents (see Figure 3, Table 5).   The blue segments in the histograms highlight the orientation, presence, and nearness of dimeric paracetamol molecules to the orthogonal synthons as represented in form I, and the red parts of histograms refers to the nearness of the orientation of paracetamol molecules to the synthons that are representative of form II. The profile of hydrogen bonds varied significantly with respect to each other. For chloroform, the percentage distribution within the range of 45-135° was 66%. The peak in the middle at 90° signifies orthogonal positions of participating molecules when hydrogen bonds were detected. Conversely, in case of paracetamol in water simulation results, the distribution of hydrogen bond showed only 14.14% hydrogen bond within a range of 45-135° and the dip in the middle indicated the absence of this orthogonal orientation and hinted towards the parallel positions found in form II. When water was used as a solvent, paracetamol crystallized in form II. Interestingly, the closely related solvents ethanol and methanol yielded different polymorphs of paracetamol when used as a solvent. The profiles of the hydrogen bonds in the simulations looked similar, although the percentage of hydrogen bonds falling within the range of 45-135° was very different, for ethanol it The blue segments in the histograms highlight the orientation, presence, and nearness of dimeric paracetamol molecules to the orthogonal synthons as represented in form I, and the red parts of histograms refers to the nearness of the orientation of paracetamol molecules to the synthons that are representative of form II. The profile of hydrogen bonds varied significantly with respect to each other. For chloroform, the percentage distribution within the range of 45-135 • was 66%. The peak in the middle at 90 • signifies orthogonal positions of participating molecules when hydrogen bonds were detected. Conversely, in case of paracetamol in water simulation results, the distribution of hydrogen bond showed only 14.14% hydrogen bond within a range of 45-135 • and the dip in the middle indicated the absence of this orthogonal orientation and hinted towards the parallel positions found in form II. When water was used as a solvent, paracetamol crystallized in form II. Interestingly, the closely related solvents ethanol and methanol yielded different polymorphs of paracetamol when used as a solvent. The profiles of the hydrogen bonds in the simulations looked similar, although the percentage of hydrogen bonds falling within the range of 45-135 • was very different, for ethanol it was 64.26%, and for methanol, it was 46.28% to still term them as form I and form II. The chemical nature, miscibility within each other, the total number of hydrogen bonds in the simulations, and profile of the histograms, and eventually different polymorphs returned by these two solvents made ethanol and methanol intriguing cases to explore further. In addition to this, paracetamol itself is an aminophenol. Therefore, we decided to dope each solvent with another to check the effect on the crystallization in the wet lab. The results are tabulated in Table 6.

MD Simulation
The different solvents in the MD trajectory were distinct in its effect on the self-association of paracetamol molecules. The simplification of larger, complex paracetamol superstructure to a synthon is a practicable methodology in analyzing the minimal unit responsible for the effect that could be the cause of a completely different end result, a different polymorph in case of paracetamol. Such a structural link suggested that at least at the dimer level of paracetamol, the nature of the self-associate and its intermolecular bonding can be an essential factor in the crystal nucleation process. An added investigative procedure to restrain/stabilize dimeric synthons and continuous evaluative aggregates formation can be next logical step in understanding role of the dimers in crystallization.

Wet Lab Experiments
Ethanol and methanol as the choice of solvent had an effect where the crystallization was driven to form either polymorph I or II. When we slowly started doping ethanol as solvent by methanol, the system became unstable. At 95% of methanol, high concentrations of methanol doping the system collapsed and no definitive polymorph could be obtained from the solvent system. This could also be interpreted as, even at very low concentration of ethanol, the ability of methanol to drive the process towards form II was highly compromised. The methanol behaved as a weak solvent in this case, since similar behavior was not observed at high concentrations of ethanol where methanol was present in a low fraction of the whole system. This strong-arming of the solvents with respect to each other was clearly observed in the case of the solvent pair of ethanol and methanol.

Conclusions
In our efforts to understand solvent effects on the crystallization of paracetamol. The molecular dynamics simulations indicated that solvent-solute interactions influence the preferred hydrogen bonding within the paracetamol molecules. The total number of hydrogen bonds and the respective positions of the contributing paracetamol molecules in a hydrogen bond was affected by choice of the solvent. Pattern recognition of synthons originating from crystallized polymorphs in the MD simulation trajectories can serve as an important tool in understanding the effects of solvents on the crystallization of drug molecules, and hence exert control on such processes.
The experiments with mixtures of solvents clearly showed a difference in solvent strengths leading crystallization process towards a polymorph. The experimental results concurred with the predictions made based on the number of hydrogen bonds in pure solvents. These types of studies can pave a path for differentiating solvents and ranking them for the desired outcome for a particular chemical class of compounds.