Recent Advances in Co-Former Screening and Formation Prediction of Multicomponent Solid Forms of Low Molecular Weight Drugs

Multicomponent solid forms of low molecular weight drugs, such as co-crystals, salts, and co-amorphous systems, are a result of the combination of an active pharmaceutical ingredient (API) with a pharmaceutically acceptable co-former. These solid forms can enhance the physicochemical and pharmacokinetic properties of APIs, making them increasingly interesting and important in recent decades. Nevertheless, predicting the formation of API multicomponent solid forms in the early stages of formulation development can be challenging, as it often requires significant time and resources. To address this, empirical and computational methods have been developed to help screen for potential co-formers more efficiently and accurately, thus reducing the number of laboratory experiments needed. This review provides a comprehensive overview of current screening and prediction methods for the formation of API multicomponent solid forms, covering both crystalline states (co-crystals and salts) and amorphous forms (co-amorphous). Furthermore, it discusses recent advances and emerging trends in prediction methods, with a particular focus on artificial intelligence.


Introduction
Drug absorption after oral administration of active pharmaceutical ingredients (APIs) inter alia depends on their physicochemical properties (e.g., polymorphic form, drug solubility), their dose, and the local environment within the gastrointestinal tract.It has been reported that up to 90% of the currently developed APIs and about 40% of the approved drugs have poor biopharmaceutical properties as a consequence of their insufficient solubility in water [1,2].The pharmaceutical industry is thus investing considerable efforts in searching for strategies to improve the solubility and bioavailability of APIs, without compromising their effectiveness.Several strategies [3], including nano-and micro-based drug delivery systems (e.g., extracellular vesicles [4], nanospheres and micelles [5], and nano-and microemulsions [6]), modified release solid dosage forms (e.g., in tablets and capsules [7]), and crystal engineering [8] (e.g., co-crystals [9] and co-amorphous forms) have been successfully investigated and partly implemented to improve the solubility and dissolution rate of APIs belonging to the biopharmaceutics classification system classes II and IV [10].
In 1989, Desiraju defined the term crystal engineering as "the comprehension of intermolecular interactions within crystal packing, and the use of this knowledge to engineer novel solid materials possessing specific physical and chemical characteristics" [11].Crystal In 1989, Desiraju defined the term crystal engineering as "the comprehension of intermolecular interactions within crystal packing, and the use of this knowledge to engineer novel solid materials possessing specific physical and chemical characteristics" [11].Crystal engineering is a powerful tool in the design of crystalline structures presenting different packing and intermolecular interactions.These crystalline structures are composed of one or more types of molecules (multicomponent solid forms).[12] Among these solid forms, co-crystals, salts, and co-amorphous systems (Figure 1) have been considered important in improving the physicochemical properties of APIs [13][14][15].In 2018, the U.S. Food and Drug Administration (FDA) defined pharmaceutical cocrystals as crystalline materials composed of a neutral API and a second neutral molecule generally defined as a co-former (Figure 1).Co-crystals are formed in a stoichiometric ratio where the molecules interact via non-covalent interactions, such as hydrogen bonds, van der Waals interactions, π•••π stacking, and halogen bonds, to form crystalline structures [13].According to the FDA, this is the only difference between a co-crystal and a salt, as, in a pharmaceutical salt, a proton transfer occurs between an ionizable API and the coformer and both molecules interact in a stoichiometric way via charge-assisted hydrogen bond interactions (Figure 1) [16].Solvates and hydrates are crystalline materials in which solvent molecules are present in the crystal lattice.Solvent molecules can thus be present in the crystalline lattice of co-crystals and salts, forming solvated/hydrated co-crystals and solvated/hydrated salts [16].
Amorphous solid dispersion technology has evolved since the 1960s, with polymers used as carriers.However, polymers have limitations, such as hygroscopicity and low drug loading.In 2009, Chieng et al. [17] introduced a term called "co-amorphous", for amorphous solid dispersions that replace polymers with low molecular weight compounds.Co-amorphous refers to a homogenous, single-phase amorphous system containing two or more low molecular weight components interacting with each other in a non-periodic way, via, e.g., charge or non-charge assisted hydrogen bond interactions [18].Based on the type of co-former used, co-amorphous systems can be divided into APIexcipient and API-API co-amorphous systems [19].In API-excipient co-amorphous systems, saccharin, nicotinamide, amino acids, dipeptides and carboxylic acids can be used as excipients.[20] On the other hand, in API-API co-amorphous systems, APIs can be selected according to their similar or complementary pharmacological properties [21].
Suitable co-formers are crucial to designing stable co-crystal, salts, and coamorphous systems with desirable properties, such as high solubility, fast dissolution rate, good diffusion permeability, high physicochemical stability and a possible synergistic pharmacological effect [22].Over the last few decades, high-throughput screening experiments have been used to predict the formation of API multicomponent solid forms.This trial and error approach is time consuming and requires a considerable consumption of materials and resources [23].In 2018, the U.S. Food and Drug Administration (FDA) defined pharmaceutical cocrystals as crystalline materials composed of a neutral API and a second neutral molecule generally defined as a co-former (Figure 1).Co-crystals are formed in a stoichiometric ratio where the molecules interact via non-covalent interactions, such as hydrogen bonds, van der Waals interactions, π•••π stacking, and halogen bonds, to form crystalline structures [13].According to the FDA, this is the only difference between a co-crystal and a salt, as, in a pharmaceutical salt, a proton transfer occurs between an ionizable API and the co-former and both molecules interact in a stoichiometric way via charge-assisted hydrogen bond interactions (Figure 1) [16].Solvates and hydrates are crystalline materials in which solvent molecules are present in the crystal lattice.Solvent molecules can thus be present in the crystalline lattice of co-crystals and salts, forming solvated/hydrated co-crystals and solvated/hydrated salts [16].
Amorphous solid dispersion technology has evolved since the 1960s, with polymers used as carriers.However, polymers have limitations, such as hygroscopicity and low drug loading.In 2009, Chieng et al. [17] introduced a term called "co-amorphous", for amorphous solid dispersions that replace polymers with low molecular weight compounds.Co-amorphous refers to a homogenous, single-phase amorphous system containing two or more low molecular weight components interacting with each other in a non-periodic way, via, e.g., charge or non-charge assisted hydrogen bond interactions [18].Based on the type of co-former used, co-amorphous systems can be divided into API-excipient and API-API co-amorphous systems [19].In API-excipient co-amorphous systems, saccharin, nicotinamide, amino acids, dipeptides and carboxylic acids can be used as excipients.[20] On the other hand, in API-API co-amorphous systems, APIs can be selected according to their similar or complementary pharmacological properties [21].
Suitable co-formers are crucial to designing stable co-crystal, salts, and co-amorphous systems with desirable properties, such as high solubility, fast dissolution rate, good diffusion permeability, high physicochemical stability and a possible synergistic pharmacological effect [22].Over the last few decades, high-throughput screening experiments have been used to predict the formation of API multicomponent solid forms.This trial and error approach is time consuming and requires a considerable consumption of materials and resources [23].

•
Quantitative analysis of MEPs and periodic density functional theory calculations were used to examine the relationship between the donor/acceptor groups in TDZH/co-former molecules and the hydrogen bond pattern in multicomponent crystal forms. [37]

•
The results met well with the ∆pKa rule: below 0, co-crystals are formed; above 3, salts are formed; and in the middle range, the formation of both is favored.

•
The PLS-DA model effectively separated co-amorphous and non-co-amorphous samples.

•
Non-polar side chain amino acids were the most effective co-formers for the formation of co-amorphous systems, while polar amino acids were the least successful.

∆pK a Rule
The ∆pK a Rule is a simple approach to predict the formation of co-crystals or salts [63,68].The difference in pK a (∆pK a ) between acidic/basic APIs and basic/acidic co-formers is calculated as follows: ∆pK a = pK a (base) − pK a (acid) ∆pK a can be used to estimate the tendency for proton transfer between a given API and co-former.When ∆pK a is higher than 3, a significant difference in the acidity or basicity of the API and co-former is obtained, indicating a preferable formation of a salt.When ∆pK a is lower than 0, the acidity or basicity of the API and co-former are similar and thus the formation of a co-crystal is expected where non-charge-assisted hydrogen bond interactions occur between the API and co-former.However, when ∆pK a is between 0 and 3, the difference in acidity or basicity is not large enough to clearly favor the formation of salt over a co-crystal or vice versa.In such cases, the ∆pK a rule is not the most reliable approach to predict the formation of a salt or co-crystal [9,69].A combination of experimental and computational methods is thus necessary to determine the proton location.These methods include nuclear magnetic resonance spectroscopy, X-ray crystallography, and computational modeling using density functional theory or molecular dynamics simulations [70].The ∆pK a rule was validated by Cruz-Cabeza using 6465 crystalline structures from the CSD [71].Multicomponent crystalline solid forms with ionized or non-ionized acid-base pairs are only observed when ∆pK a is greater than 4 or less than -1, respectively.When ∆pK a is between -1 and 4, ∆pK a and the likelihood of proton transfer showed a linear relationship (Figure 3).A "salt-co-crystal continuum" exists when the proton position is ambiguous (∆pK a between 0 and 3), which is between the two extremes of salt and co-crystal [63].In such cases, the ∆pK a rule is not applicable for predicting the formation of salts or co-crystals.Childs et al. [63], investigated the propensity for forming co-crystals or salts with multicomponent pairs with ∆pK a lower than 3.They prepared a total of 20 multicomponent solid forms of theophylline and guest molecules, which included 13 co-crystals, five salts, and two within the salt-cocrystal continuum.It is important to note that the formation of salts or co-crystals is influenced by several factors, including the solvent used (and its polarity), temperature, API or coformer concentration, and crystal packing interactions.Therefore, the values of ∆pK a alone cannot predict with certainty the formation of salts or co-crystals, but they can provide useful information about the likelihood of their formation.

Supramolecular Synthons
Recognizing supramolecular interactions is crucial for designing crystals.After introducing the term "crystal engineering" in 1989, Desiraju [72] introduced the term "supramolecular synthons" and defined them as structural entities within supermolecules (i.e., complexes of two or more molecules that are non-covalently bonded) that can be created and/or arranged through intermolecular interactions using existing or feasible

Supramolecular Synthons
Recognizing supramolecular interactions is crucial for designing crystals.After introducing the term "crystal engineering" in 1989, Desiraju [72] introduced the term "supramolecular synthons" and defined them as structural entities within supermolecules (i.e., complexes of two or more molecules that are non-covalently bonded) that can be created and/or arranged through intermolecular interactions using existing or feasible synthetic methods.Supramolecular synthons can serve as a design strategy for controlling the self-assembly of molecules through non-covalent interactions in the solid state.These synthons play a fundamental role in the formation of co-crystals because they act as building blocks in supramolecular chemistry, guiding the arrangement and organization of molecules in the solid state.Supramolecular synthons can be divided into two categories: homosynthons and heterosynthons (Figure 4).Supramolecular homosynthons consist of the same complementary functional groups that can form self-association motifs, such as dimers or chains, including acid-acid [73] and amide-amide [74] interactions.Supramolecular heterosynthons consist of diverse but complementary functional groups, including acid-amide [75], acid-hydroxyl [76], hydroxyl-pyridine [77], acid-pyridine [28], acid-Noxide [70], amide-N-oxide [78], sulfonamide-N-oxide [79], and sulfonamide-amide [80] interactions.In a study conducted by the CSD, the hierarchy of supramolecular heterosynthons, involving carboxylic acids and alcohols, was evaluated in competitive environments, i.e., considering that both carboxylic acid and hydroxyl groups compete to form hydrogen bond motifs [28] were chosen as co-former candidates because hydroxyl and carbonyl groups of carboxylic acids can act as hydrogen bond donor and acceptor, respectively.Three dicarboxylic acids including malonic acid, glutaric acid and pimelic acid were found to form co-crystals with REG.Based on the knowledge that nicotinamide and isonicotinamide commonly form heterosynthons with carboxylic acids, Das et al. [32] used nitrogen-containing carboxylic acids, including 3,5-pyrazole dicaboxylic acid, dipicolinic acid, or quinolinic acid as coformers and successfully prepared new multicomponent crystal forms.Four different supramolecular heterosynthons were found in nicotinamide and isonicotinamide cocrystals with 3,5-pyrazole dicarboxylic acid (Figure 5).In co-former screening for minoxidil, Deng et al.The use of supramolecular synthons has been explored only in multicomponent crystal forms [31,32,70].Since the supramolecular synthon principle is simple and does not require any calculations, it has been most frequently used as a tool to guide co-former selection.For example, when screening co-crystals of regorafenib (REG), it was found that REG has several hydrogen bond sites including hydrogen bond donors (amino groups), and hydrogen bond acceptors (carbonyl and pyridine groups) [31].Dicarboxylic acids were chosen as co-former candidates because hydroxyl and carbonyl groups of carboxylic acids can act as hydrogen bond donor and acceptor, respectively.Three dicarboxylic acids including malonic acid, glutaric acid and pimelic acid were found to form co-crystals with REG.Based on the knowledge that nicotinamide and isonicotinamide commonly form heterosynthons with carboxylic acids, Das et al. [32] used nitrogen-containing carboxylic acids, including 3,5-pyrazole dicaboxylic acid, dipicolinic acid, or quinolinic acid as co-formers and successfully prepared new multicomponent crystal forms.Four different supramolecular heterosynthons were found in nicotinamide and isonicotinamide co-crystals with 3,5-pyrazole dicarboxylic acid (Figure 5).In co-former screening for minoxidil, Deng et al. [70]

Molecular Electrostatic Potential Surfaces (MEPs)
The calculated MEPs generated using the density functional theory in the gas phase are used to identify surface site interaction points and predict electrostatic interactions at the surface of the molecules [81].The strength of hydrogen bond donors and acceptors can be ranked according to the MEPs [82], which have also been further used in the design of ternary multicomponent crystal forms [83].The hydrogen bond donor (α) and acceptor (β) are generated from the maxima and minima of the MEP and calculated using Equations ( 2) and (3), respectively.0.0000162MEP 0.00962MEP where MEPmin and MEPmax are the local minima and maxima on the MEPs (unit: kJ•mol −1 ).The prediction of multicomponent crystal forms based on MEPs considers that all hydrogen bond sites on the surface of a molecule are not restricted by the internal molecular structure and are free to interact with other molecules in the solid-state

Molecular Electrostatic Potential Surfaces (MEPs)
The calculated MEPs generated using the density functional theory in the gas phase are used to identify surface site interaction points and predict electrostatic interactions at the surface of the molecules [81].The strength of hydrogen bond donors and acceptors can be ranked according to the MEPs [82], which have also been further used in the design of ternary multicomponent crystal forms [83].The hydrogen bond donor (α) and acceptor (β) are generated from the maxima and minima of the MEP and calculated using Equations ( 2) and (3), respectively.
where MEP min and MEP max are the local minima and maxima on the MEPs (unit: kJ•mol −1 ).The prediction of multicomponent crystal forms based on MEPs considers that all hydrogen bond sites on the surface of a molecule are not restricted by the internal molecular structure and are free to interact with other molecules in the solid-state environment.The molecular arrangement, steric constraints and packing effects are not taken into account [81].The interaction site pairing energy E is the sum of all contacts across the surface of every molecule in the crystal, which can be calculated using Equation ( 4) [84].
where α i are hydrogen bond donor sites, and β j are hydrogen bond acceptor sites.Positive or negative unpaired sites locate low electrostatic potential gaps or regions, making no contribution to the overall electrostatic interaction energy [81].
According to Etter, the first most positive α i interacts with the first most negative β j , the second most positive α i with the second most negative β j , and so forth until all of the interaction sites are covered [85].Based on this theory, Musumeci et al. [81] proposed that the probability of forming a multicomponent crystal form can be predicted by the difference in the interaction site pairing energy (∆E) between the total E of pure components and the E of the multicomponent crystal forms, using the following Equation (5).
where E crys , E 1 and E 2 are the interaction site pairing energies of a multicomponent crystal of stoichiometry 1 n 2 m , and the pure component solids, 1 and 2, respectively.The interaction site pairing energy of a multicomponent crystal is calculated in the same way as for a pure component solid.A high value of ∆E indicates a stronger interaction between two different components and a higher probability of forming a multicomponent crystal form.The minimum value of ∆E is 0 kJ•mol −1 , indicating that the formation of multicomponent crystal forms must increase the interaction energy [81].The specific cut-off value depends on the system involved.
The prediction reliability of MEPs has until now only been demonstrated for multicomponent crystal forms.For example, Pagliari et al. [86] applied MEP plots to illustrate the complementarity between amide sites and aromatic rings, which clarified how co-crystals are formed.Musumeci et al. [81] used caffeine and carbamazepine as model APIs and about 1000 co-former candidates to examine the reliability of the prediction method.The results demonstrated that when ∆E was over 11 kJ•mol −1 , the possibility of co-crystal formation is larger than 50%.Grecu et al. [84] successfully validated MEPs as a virtual co-crystal screening tool using reported cases of 18 APIs.The co-former candidates were ranked according to their ∆E values and the ones presenting a larger value of ∆E were consistent with the experimental results in most cases.Furthermore, the MEPs method showed little difference when compared to the COSMO-based methods (see below).The MEPs method can also be an effective tool to (i) investigate the driving force of co-crystals/salts formation and explain the ratio of APIs and co-formers in multicomponent crystal forms [65,70], and (ii) elucidate the topology of hydrogen bonds and identify the factors contributing to any observed disorder in a crystal lattice [37].Apart from hydrogen bond interaction, the MEPs method was also explored to give insight into halogen bond [87] and chalcogen bond [88] formation, which belongs to σ hole interactions and is coming into more focus recently.For instance, in co-crystals of 1,4-diiodotetrafluorobenzene and the isomeric n-pyridinealdazines (n = 2, 3 and 4), the I•••N halogen bond interaction is the primary interaction [87].The σand π-holes of the MEPs in 1,4-diiodotetrafluorobenzene showed high potential and were used as halogen bond donors.The halogen bond strength order for n-pyridinealdazines was 4 > 3 > 2 according to the V s,min (minimum potentials) values located near the pyridine-N atom for n = 2, 3 and 4 (Figure 6).Wzgarda-Raj et al. [88] synthesized four new multicomponent crystals of trithiocyanuric acid with pyridine N-oxide derivatives, where N-H•••S hydrogen bonds were observed to form R 2 2 (8) synthons.Overall, the validated approach provides a promising framework for virtual multicomponent crystal forms screening and could potentially be applied in drug discovery.

Hydrogen Bond Propensity (HBP)
HBP was first developed by Galek et al. [89] as a logit (i.e., probabilistic) hydrogen bond propensity method to quantify the probability of hydrogen bond formation and used in the experimental polymorphic screening of ritonavir.This method was further used to successfully screen for multicomponent crystal forms of lamotrigine [90].Based on the HBP model, it is important to recognize every potential hydrogen bond by identifying the donors and acceptors, which are responsible for the formation of multicomponent crystal forms.
HBP was developed by the Cambridge Crystallographic Data Centre and has been integrated into the Mercury program [91] to identify donors and acceptors forming usual or unusual hydrogen bonds.Molecules from the CSD with certain functional groups were used to prepare a dataset to train statistical models [41].When predicting cocrystallization between an API and a co-former, both homomeric and heteromeric interactions are taken into account.If HBPAPI-co-former is larger than HBPAPI-API and HBPcoformer-co-former, co-crystallization is likely to occur [92].In Figure 7, A and B are two different molecules.A-A and B-B are homomeric interactions within A and B molecules, respectively.A-B and B-A are two different heteromeric interactions between A and B molecules.The Δpropensity value is calculated as the difference of largest heteromeric interactions and largest homomeric interactions as shown in Equation (6).

Hydrogen Bond Propensity (HBP)
HBP was first developed by Galek et al. [89] as a logit (i.e., probabilistic) hydrogen bond propensity method to quantify the probability of hydrogen bond formation and used in the experimental polymorphic screening of ritonavir.This method was further used to successfully screen for multicomponent crystal forms of lamotrigine [90].Based on the HBP model, it is important to recognize every potential hydrogen bond by identifying the donors and acceptors, which are responsible for the formation of multicomponent crystal forms.
HBP was developed by the Cambridge Crystallographic Data Centre and has been integrated into the Mercury program [91] to identify donors and acceptors forming usual or unusual hydrogen bonds.Molecules from the CSD with certain functional groups were used to prepare a dataset to train statistical models [41].When predicting co-crystallization between an API and a co-former, both homomeric and heteromeric interactions are taken into account.If HBP API-co-former is larger than HBP API-API and HBP co-former-co-former , cocrystallization is likely to occur [92].In Figure 7, A and B are two different molecules.A-A and B-B are homomeric interactions within A and B molecules, respectively.A-B and B-A are two different heteromeric interactions between A and B molecules.The ∆ propensity value is calculated as the difference of largest heteromeric interactions and largest homomeric interactions as shown in Equation (6).
Figure 7. HBP approach for screening of multicomponent crystal forms.
The HBP method gives information about the possibility of forming a specific hydrogen bond, and it depends on how frequently the interaction occurs with respect to other interactions in a given fitting data set.In the HBP calculation for the system lenalidomide-nicotinamide [42], the -NH2 of nicotinamide group and the carbonyl group of lenalidomide were found to have a high propensity for heteromeric interactions (0.94), while homomeric interactions between lenalidomide-lenalidomide have a propensity of 0.84 for -NH-with carbonyl group, and 0.82 between nicotinamide-nicotinamide for -NH2 with carbonyl group (Figure 8).Heteromeric interactions were preferred as indicated by a Δpropensity value of 0.1 and the co-crystal was successfully prepared experimentally.Majumde et al. [41] applied HBP in predicting indomethacin-nicotinamide (1:1) cocrystals, finding that their hydrogen bond motifs (N3-H are the most likely donor-acceptor combinations predicted by HBP.Sarkar et al. [92] used HBP to predict the co-crystal formation between six model APIs and 25 possible coformers.A success rate of 92-95% was obtained, indicating an excellent hydrogen bond affinity between APIs and co-formers.In another study, Sarkar et al. [93] used HBP, molecular complementarity, and hydrogen bond energy to predict co-crystal formation between six APIs and 42 potential co-formers.The combination of HBP and molecular complementarity allowed the achievement of an overall accuracy of 81%.Since HBP is a hydrogen bond-based method, its application is limited to multicomponent solid forms where hydrogen bonds are the dominant interaction [94].

Lattice Energy
Lattice energy is another tool to predict multicomponent crystal forms.The energy difference between multicomponent crystal forms and pure individual components can be used as an indication of whether co-crystallization is expected to occur spontaneously, The HBP method gives information about the possibility of forming a specific hydrogen bond, and it depends on how frequently the interaction occurs with respect to other interactions in a given fitting data set.In the HBP calculation for the system lenalidomidenicotinamide [42], the -NH 2 of nicotinamide group and the carbonyl group of lenalidomide were found to have a high propensity for heteromeric interactions (0.94), while homomeric interactions between lenalidomide-lenalidomide have a propensity of 0.84 for -NHwith carbonyl group, and 0.82 between nicotinamide-nicotinamide for -NH 2 with carbonyl group (Figure 8).Heteromeric interactions were preferred as indicated by a ∆ propensity value of 0.1 and the co-crystal was successfully prepared experimentally.Majumde et al. [41] applied HBP in predicting indomethacin-nicotinamide (1:1) co-crystals, finding that their hydrogen bond motifs (N3-H•••O5, N3-H•••O4 and O3-H•••N2), are the most likely donoracceptor combinations predicted by HBP.Sarkar et al. [92] used HBP to predict the co-crystal formation between six model APIs and 25 possible co-formers.A success rate of 92-95% was obtained, indicating an excellent hydrogen bond affinity between APIs and co-formers.In another study, Sarkar et al. [93] used HBP, molecular complementarity, and hydrogen bond energy to predict co-crystal formation between six APIs and 42 potential co-formers.The combination of HBP and molecular complementarity allowed the achievement of an overall accuracy of 81%.Since HBP is a hydrogen bond-based method, its application is limited to multicomponent solid forms where hydrogen bonds are the dominant interaction [94].The HBP method gives information about the possibility of forming a specific hydrogen bond, and it depends on how frequently the interaction occurs with respect to other interactions in a given fitting data set.In the HBP calculation for the system lenalidomide-nicotinamide [42], the -NH2 of nicotinamide group and the carbonyl group of lenalidomide were found to have a high propensity for heteromeric interactions (0.94), while homomeric interactions between lenalidomide-lenalidomide have a propensity of 0.84 for -NH-with carbonyl group, and 0.82 between nicotinamide-nicotinamide for -NH2 with carbonyl group (Figure 8).Heteromeric interactions were preferred as indicated by a Δpropensity value of 0.1 and the co-crystal was successfully prepared experimentally.Majumde et al. [41] applied HBP in predicting indomethacin-nicotinamide (1:1) cocrystals, finding that their hydrogen bond motifs (N3-H are the most likely donor-acceptor combinations predicted by HBP.Sarkar et al. [92] used HBP to predict the co-crystal formation between six model APIs and 25 possible coformers.A success rate of 92-95% was obtained, indicating an excellent hydrogen bond affinity between APIs and co-formers.In another study, Sarkar et al. [93] used HBP, molecular complementarity, and hydrogen bond energy to predict co-crystal formation between six APIs and 42 potential co-formers.The combination of HBP and molecular complementarity allowed the achievement of an overall accuracy of 81%.Since HBP is a hydrogen bond-based method, its application is limited to multicomponent solid forms where hydrogen bonds are the dominant interaction [94].

Lattice Energy
Lattice energy is another tool to predict multicomponent crystal forms.The energy difference between multicomponent crystal forms and pure individual components can be used as an indication of whether co-crystallization is expected to occur spontaneously,

Lattice Energy
Lattice energy is another tool to predict multicomponent crystal forms.The energy difference between multicomponent crystal forms and pure individual components can be used as an indication of whether co-crystallization is expected to occur spontaneously, i.e., if it is thermodynamically favored compared to the crystallization of each starting material separately [95].Considering that co-crystallization is commonly found to be a thermodynamically favorable process, its lattice energy is higher compared to the corresponding individual components [69,95,96].The difference in lattice energies ∆E latt can be calculated as follows: where E latt (A m B n ) is the lattice energy of multicomponent crystal forms A m B n consisting of molecules A and B in a stoichiometric ratio m:n, and E latt (A) and E latt (B) are the lattice energies of the pure components A and B, respectively.The probability of forming a multicomponent crystal form is higher when the ∆E latt is more negative.It is unlikely that multicomponent crystal forms are formed if ∆E latt has a positive value.Distinct from other methods, this prediction technique does not make any assumptions about hydrogen bonds and relies solely on the lattice energy.Factors such as temperature, pressure, solvent effects, and kinetics are not considered in the calculations [96].
The E latt of the crystal can be calculated by three methods [97]: the first one is using a set of programs for crystallography, FlexCryst [98], to calculate free energy G.The second approach is using cohesive energy [99], and the last one is using total energy resulting from noncovalent, pairwise interactions between the molecule and its surrounding molecules [100].Chan et al. [96] performed lattice energy calculations and investigated the thermodynamic stability of 102 multicomponent crystal forms containing nicotinamide, isonicotinamide, and picolinamide, 99 of which (more than 97%) were consistent with the observed tendency of the compound to crystallize.Kuleshova and co-workers [101] used the FlexCryst program suite to calculate the free energy and determine the relative stability of co-crystals of flavonoids and their pure crystal forms taken from the CSD.It was found that the lattice energy calculation was a valuable tool for in silico screening of co-crystal formation and stability, which can further be used to estimate their relative solubility.Sun and co-workers [102] proposed virtual co-former screening approaches, taking into account the lattice packing contributions of crystals to screen co-formers for indomethacin and paracetamol.If the lattice energy difference is ranked from smallest to largest, successful co-crystal formation can be found in the first six values of the lists (Figure 9).Lattice energy can also be used to determine the stable form of co-crystals.Surov et al. [103] compared the theoretical lattice energies and found that in multicomponent crystal forms of fluconazole with 4-hydroxybenzoic acid, hydrated crystal forms are more energetically favorable than the anhydrous co-crystals.Moreover, based on energy calculations, Vener and co-workers [97] found that the range of supramolecular synthons is approximately ~80 to ~30 kJ/mol, with a decreasing order of strength as follows: acid-amide > acid-pyridine > hydroxyl-acid > amide-amide > hydroxyl-pyridine.

Molecular Complementarity (MC, Fábián's Method)
Molecular complementarity (MC) was first introduced by Fábián in 2009 [104] to investigate the molecular characteristics that affect the formation of multicomponent crystal forms.In his study, a statistical analysis was carried out using 131 descriptors of 1949 molecules.It was pointed out that molecules with similar properties, especially in molecular shape and polarity, tend to have a higher propensity to form stable multicomponent crystal forms.Five numerical descriptors including three shape descriptors (S-axis, S/L axis, and M/L axis) and two polarity descriptors (fraction of nitrogen and oxygen atoms, and dipole moment) were considered to be important in the formation of multicomponent crystal forms.S-axis is the length of the short axis, S/L axis is the short/long axis ratio, and M/L axis is the medium/long axis ratio, while L, M and S refer to the three unequal dimensions of a rectangular crystal cell, as shown in Figure 10 [93].The MC method has been developed by the Cambridge Crystallographic Data Centre and incorporated into the Mercury software (version 3.7) [94].Every descriptor has a criterion to indicate "PASS" or "FAIL".A "PASS" indicates that a multicomponent crystal form was successfully formed, while a "FAIL" indicates that no multicomponent crystal forms were formed [104,105].A multicomponent crystal form is likely to be formed only when all five descriptors demonstrate a "PASS."

Molecular Complementarity (MC, Fábián's Method)
Molecular complementarity (MC) was first introduced by Fábián in 2009 [104] to investigate the molecular characteristics that affect the formation of multicomponent crystal forms.In his study, a statistical analysis was carried out using 131 descriptors of 1949 molecules.It was pointed out that molecules with similar properties, especially in molecular shape and polarity, tend to have a higher propensity to form stable multicomponent crystal forms.Five numerical descriptors including three shape descriptors (S-axis, S/L axis, and M/L axis) and two polarity descriptors (fraction of nitrogen and oxygen atoms, and dipole moment) were considered to be important in the formation of multicomponent crystal forms.S-axis is the length of the short axis, S/L axis is the short/long axis ratio, and M/L axis is the medium/long axis ratio, while L, M and S refer to the three unequal dimensions of a rectangular crystal cell, as shown in Figure 10 [93].The MC method has been developed by the Cambridge Crystallographic Data Centre and incorporated into the Mercury software (version 3.7) [94].Every descriptor has a criterion to indicate "PASS" or "FAIL".A "PASS" indicates that a multicomponent crystal form was successfully formed, while a "FAIL" indicates that no multicomponent crystal forms were formed [104,105].A multicomponent crystal form is likely to be formed only when all five descriptors demonstrate a "PASS." MC is mostly used as a preliminary screening tool instead of being used solely to determine the formation of multicomponent crystal forms.For example, Li et al. [106] used both COSMO-RS and MC in screening the formation of new multicomponent crystal forms of 2,4-dichlorophenoxyacetic acid.In total, 25 out of the 53 co-former candidates were identified to be potential co-formers, among which 20 co-formers were experimentally successful in forming new multicomponent crystal forms with 2,4dichlorophenoxyacetic acid (Figure 11).Wu et al. [65] evaluated COSMO-RS, MC, HSP and their combinations in screening multicomponent crystal forms of 2-amino-4,6dimethoxypyrimidine with 63 components.The overall successful rate of MC was 69.8%, and the best outcomes were obtained when using MC and COSMO-RS, with an overall

Medium axis
Short axis (S) MC is mostly used as a preliminary screening tool instead of being used solely to determine the formation of multicomponent crystal forms.For example, Li et al. [106] used both COSMO-RS and MC in screening the formation of new multicomponent crystal forms of 2,4-dichlorophenoxyacetic acid.In total, 25 out of the 53 co-former candidates were identified to be potential co-formers, among which 20 co-formers were experimentally successful in forming new multicomponent crystal forms with 2,4-dichlorophenoxyacetic acid (Figure 11).Wu et al. [65] evaluated COSMO-RS, MC, HSP and their combinations in screening multicomponent crystal forms of 2-amino-4,6-dimethoxypyrimidine with 63 components.The overall successful rate of MC was 69.8%, and the best outcomes were obtained when using MC and COSMO-RS, with an overall success rate of 85.7%.Wu et al. [64] also used a combined method (COSMO-RS and MC) to screen multicomponent crystal forms for a pesticide, pymetrozine, with 39 co-former candidates.The calculation resulted in 13 promising co-formers, and it was discovered that nine of them produced novel solid phases.

Hansen Solubility Parameter (HSP)
The solubility parameter theory was proposed by Hildebrand and Scott in 1950 [107].According to this theory, the cohesive energy indicates the total interactions in a material, including hydrogen bonds, van der Waals forces, covalent bonds and ionic bonds.Cohesive energy per unit volume, i.e., the cohesive energy density (CED), is an essential Medium axis (M) Short axis (S)

Hansen Solubility Parameter (HSP)
The solubility parameter theory was proposed by Hildebrand and Scott in 1950 [107].According to this theory, the cohesive energy indicates the total interactions in a material, including hydrogen bonds, van der Waals forces, covalent bonds and ionic bonds.Cohesive energy per unit volume, i.e., the cohesive energy density (CED), is an essential parameter used in pharmaceutical research to predict the physical and chemical properties of drugs, excipients and carriers (e.g., miscibility of a drug with excipients and carriers in solid dispersions).The relationship between the CED and the solubility parameter (δ, unit: MPa 0.5 ) is shown in the following equation: where ∆E v is the energy of vaporization, and V m is the molar volume.The mutual solubility of two components is determined by their closeness in δ values.Desai and Patravale [108] applied the Hildebrand solubility parameter as one of the useful molecular descriptors to select co-formers for the formation of co-crystals with curcumin.The Hildebrand and Scott approach is based on regular solution theory and works best for non-polar molecules interacting via weak dispersion forces [109].In order to extend the application of this theory to more polar and strongly interacting systems, such as APIs, the HSP approach was developed by Hansen and divides the total solubility parameter (δ t ), or Hildebrand solubility parameter, into three partial solubility parameters, as follows [110]: where δ d are dispersion forces (atomic dispersion), δ p are 'polar' interactions, and δ h are hydrogen bonds.These partial solubility parameters are most commonly calculated based on the Hoftyzer-Van Krevelen and Fedors group contribution methods, as follows: [111] where i is the structural group within the molecule, F di , F pi , F hi and V i are the group contributions to the dispersion forces, polar forces, hydrogen bonding energy and molar volume, respectively.In 1999, Greenhalgh et al. [112] proposed employing the difference of the total solubility parameter (∆δ t ) (Equation ( 11)) to predict the miscibility between the drug and carriers, and Mohammad et al. [48] proposed a cut-off value of ∆δ t < 7 MPa 0.5 .
According to Bagley et al. [113], δ d and δ p have thermodynamically similar effects, while δ h has a different effect in nature [114], because they are interactions between hydrogen atoms and electronegative atoms such as O, N, and F. So, δ d and δ p parameters were combined as a volume-dependent solubility parameter, and δ v , and R a(v) factor were introduced as follows: The two-dimensional plot of δ v against δ h , i.e., the Bagley diagram, has been applied in different aspects, such as predicting the miscibility of two materials and the duration of drug intestinal absorption [115,116].Albers et al. [117] found that for API and polymer systems, miscibility can be achieved if R a(v) is not higher than 5.6 Mpa 0.5 .
Van Krevelen and Hoftyzer [118] introduced a three-dimensional approach to measure the difference in solubility parameters using Equation (14).
Three parameters representing hydrogen bonding, polarity, and dispersion forces can be considered as three-dimensional coordinates of HSP space points.The two components are miscible if ∆δ ≤ 5 MPa 0.5 .For the purpose of convenient visualization of the spherical, instead of using the ellipsoidal solubility plot, a modified radius (R a ) equation (Equation ( 15)) was proposed [119] as a representation of the Euclidean distance between the centrum HSP 1 (δ d1 , δ p1 , δ h1 ) and another point HSP 2 (δ d2 , δ p2 , δ h2 ) in the Hansen space (Figure 12).
( ) ( ) ( ) The HSP method was first used in co-crystal formation prediction in 2011 [48], and later, in 2020, it was used to indicate the good miscibility of tadalafil and repaglinide in co-amorphous systems [50].The reliability of the HSP as a formation prediction tool for both co-crystal/salts and co-amorphous forms has been shown in a number of articles.For example, the HSP method was used to predict the formation of multicomponent crystal forms of minoxidil (MIN) [70] with 18 co-formers, and of 2-amino-4,6dimethoxypyrimidine (MOP) [65] with 63 co-formers using the Ra(v), δ Δ and Δδt criteria.
Since there were large deviations from the experimental results for Ra(v) and δ Δ criteria, only the Δδt criterion was chosen in further discussion in both MIN and MOP cases.The overall success rates of co-former prediction by the Δδt criterion were 65.5% and 62.5% for MIN and MOP, respectively.In co-amorphous systems prediction, the Δδt value between florfenicol and oxymatrine was determined to be 3.87 MPa 0.5 , indicating good miscibility between the API and co-former, contributing to the successful formation of a coamorphous system [66].The HSP method with the van Krevelen criterion Δ δ was successful in predicting four out of six norfloxacin co-amorphous systems.The results showed better accuracy than the Greenhalgh criterion Δδt, as it could offer more accurate information about the forces present between the molecules (e.g., polar interactions, hydrogen bonding and dispersion interactions) [67].The HSP method was first used in co-crystal formation prediction in 2011 [48], and later, in 2020, it was used to indicate the good miscibility of tadalafil and repaglinide in co-amorphous systems [50].The reliability of the HSP as a formation prediction tool for both co-crystal/salts and co-amorphous forms has been shown in a number of articles.For example, the HSP method was used to predict the formation of multicomponent crystal forms of minoxidil (MIN) [70] with 18 co-formers, and of 2-amino-4,6-dimethoxypyrimidine (MOP) [65] with 63 co-formers using the R a(v) , ∆δ and ∆δ t criteria.Since there were large deviations from the experimental results for R a(v) and ∆δ criteria, only the ∆δ t criterion was chosen in further discussion in both MIN and MOP cases.The overall success rates of co-former prediction by the ∆δ t criterion were 65.5% and 62.5% for MIN and MOP, respectively.In co-amorphous systems prediction, the ∆δ t value between florfenicol and oxymatrine was determined to be 3.87 MPa 0.5 , indicating good miscibility between the API and co-former, contributing to the successful formation of a co-amorphous system [66].The HSP method with the van Krevelen criterion ∆δ was successful in predicting four out of six norfloxacin co-amorphous systems.The results showed better accuracy than the Greenhalgh criterion ∆δ t , as it could offer more accurate information about the forces present between the molecules (e.g., polar interactions, hydrogen bonding and dispersion interactions) [67].

Conductor-like Screening Model for Real Solvents (COSMO-RS)
The COSMO-RS theory was developed by Klamt and co-workers [121] to predict the thermodynamic equilibria of pure components and liquid mixtures using static thermodynamic methods based on quantum chemical calculations.According to COSMO, the solute molecule can be represented as a collection of partial charges and the surrounding solvent can be approximated as a dielectric continuum of permittivity [122].On the basis of COSMO, combined with statistical mechanics methods, Klamt [121,121] developed COSMO-RS to overcome some shortcomings of dielectric continuum solvation models.The COSMO-RS theory takes into account the molecular interactions (i.e., electrostatics, hydrogen bonding and Van der Waals interactions [123]) in calculating the chemical potentials of pure components and liquid mixtures.
COSMO-RS theory has been used in several research fields, including the prediction of solubility [122] and pKa [124], identification of suitable solvents [125], solvate or co-crystal formers [122], and the calculation of partitioning coefficients [126].It is also a possible method for screening of multicomponent crystal forms [65,70].Even though COSMO-RS theory has the capability of theoretically predicting the formation of a co-crystal or a co-amorphous system [51], the reliability of this method as a formation prediction tool only has been recently shown for multicomponent crystal forms [64,70].On the basis of the assumption that the interactions between the API and co-former in multicomponent crystal forms are comparable to those in a virtual supercooled liquid, the formation of multicomponent crystal forms can be explained by liquid phase thermodynamics without considering the long-range packing order [69].Therefore, the interaction strength of the API and co-former in the supercooled liquid state is able to predict co-crystallization, which is estimated using the excess enthalpy ∆H ex (Equation ( 16)) of the API and co-former with a given stoichiometry as compared with the pure materials [127].
where H AB is the enthalpy of the supercooled stoichiometric mixture of components A and B. H A and H B are the enthalpy of the pure components A and B, respectively, and x is the mole fraction of each component.The probability of forming a multicomponent crystal form is higher when the ∆H ex is more negative.
Although the COSMO-RS method ignores the solid-state order, it could successfully predict the formation of 21 new multicomponent crystal forms of 2-amino-4,6dimethoxypyrimidine [65] with 63 co-former candidates with an overall success rate of 84.1% when ∆H ex < −1 kcal/mol.Surprisingly, the COSMO-RS method showed up to 100% overall success rate in predicting the formation of multicomponent crystal forms of pymetrozine [64] and minoxidil [70] when ∆H ex < −3.0 kcal•mol −1 and ∆H ex < −2.00 kcal•mol −1 , respectively.Alhadid et al. [128] calculated a solid-liquid phase diagram of 5 l-menthol/xylenol eutectic systems, where menthol/3,4-xylenol at 1:2 and 2:1 ratios, and l-menthol/3,5-xylenol at 1:1 ratio could form co-crystals.The authors further investigated solid-liquid equilibria data of the l-menthol/phenol eutectic system where two co-crystals were formed in 1:2 and 2:1 ratios, and successfully used the non-random twoliquid model (an activity coefficient model frequently applied to calculate phase equilibria) and COSMO-RS models to obtain the phase diagrams [129].When exploring the possibility of forming new posaconazole co-crystals, COSMOquick selected 28 candidates from a list of approximately 140 potential co-formers.There were 13 new posaconazole multicomponent crystal forms (seven anhydrous, five hydrates, and one solvate) obtained [53].The COSMO-RS method could also detect the transferability of co-crystallization between analogous molecules.Przybyłek et al. [52] used mixing enthalpy to evaluate the co-crystal formation of two methylxanthine derivatives, theophylline and caffeine, as shown in Figure 13.The mixing enthalpy of theophylline and caffeine showed a linear relationship and phenolic acids were found to be the most promising co-formers.

Artificial Intelligence (AI)
In recent years, AI methods, especially machine learning (ML), have emerged as promising and effective strategies to predict both multicomponent crystal forms and coamorphous systems [130].This method is mainly based on molecular descriptors representing the physicochemical properties of a molecule which are computed according to its chemical structure.Figure 14 displays the ML approaches commonly used.Building data sets based on experimental data is the first step in ML [131].Molecules are characterized using artificial molecular features (such as molecular fingerprints and descriptors) or automatically extracted features (e.g., convolutional neural networks).ML algorithms are used to learn how a specific system is formed and include the creation and evaluation of models.The combination of AI with complementary experimental studies allows the leverage of the computational capabilities of AI to analyze vast amounts of data, identify patterns, and generate predictive models.Moreover, when AI predictions are complemented by experimental validations, it is possible to validate and refine these models, leading to a more comprehensive understanding of the underlying mechanisms of the formation of multicomponent solid forms.This may significantly reduce the time and cost of developing effective drugs [131], by creating a cluster analysis for co-former screening [132] and predicting the formation of multicomponent solid forms based on principal component analysis [133] and multivariate adaptive regression splines [134,135].

Artificial Intelligence (AI)
In recent years, AI methods, especially machine learning (ML), have emerged as promising and effective strategies to predict both multicomponent crystal forms and coamorphous systems [130].This method is mainly based on molecular descriptors representing the physicochemical properties of a molecule which are computed according to its chemical structure.Figure 14 displays the ML approaches commonly used.Building data sets based on experimental data is the first step in ML [131].Molecules are characterized using artificial molecular features (such as molecular fingerprints and descriptors) or automatically extracted features (e.g., convolutional neural networks).ML algorithms are used to learn how a specific system is formed and include the creation and evaluation of models.The combination of AI with complementary experimental studies allows the leverage of the computational capabilities of AI to analyze vast amounts of data, identify patterns, and generate predictive models.Moreover, when AI predictions are complemented by experimental validations, it is possible to validate and refine these models, leading to a more comprehensive understanding of the underlying mechanisms of the formation of multicomponent solid forms.This may significantly reduce the time and cost of developing effective drugs [131], by creating a cluster analysis for co-former screening [132] and predicting the formation of multicomponent solid forms based on principal component analysis [133] and multivariate adaptive regression splines [134,135].The main limitation in applying deep learning (DL) approaches in crystal engineering lies in the critical bottleneck of data availability, while the model performance also heavily relies on the key factor of data quality.The CSD comprises an extensive repository of highquality positive samples (i.e., examples of successful multicomponent solid state form formation), providing valuable support for DL applications [130].In data collection and augmentation, ML methods retrieve positive sample data from the CSD.Negative samples, however, are collected from experiments or experimental reports [130], or by computational methods such as network-based link prediction [135] and molecular The main limitation in applying deep learning (DL) approaches in crystal engineering lies in the critical bottleneck of data availability, while the model performance also heavily relies on the key factor of data quality.The CSD comprises an extensive repository of high-quality positive samples (i.e., examples of successful multicomponent solid state form formation), providing valuable support for DL applications [130].In data collection and augmentation, ML methods retrieve positive sample data from the CSD.Negative samples, however, are collected from experiments or experimental reports [130], or by computational methods such as network-based link prediction [135] and molecular similarity-based methods [57].These data can be used to train ML models.Gröls et al. [56] developed different ML algorithms based on 1000 co-crystallization events and 2083 chemical descriptors to accurately predict mechanochemical co-crystallizations.The eXtreme Gradient Boosting model was able to identify three new co-crystals of diclofenac through mechanochemistry.Similarly, they developed and compared different machine learning and deep learning algorithms based on 418 experimental amorphization cases and 2066 molecular descriptors to predict mechanochemical amorphization [59].The gradient boost model showed an accuracy of over 73% and was able to identify six novel co-amorphous systems with folic acid.Jiang et al. [130] used 6819 positive and 1052 negative samples to develop a graph neural network based on a deep learning model, as shown in Figure 15.The effectiveness of the graph neural network was shown by seven competitive models and three different and challenging out-of-sample tests, achieving prediction accuracy of higher than 96%.Wang et al. [57] developed a co-crystal prediction model based on the random forest approach.The positive samples were from the CSD and the negative ones were created by randomly combining different molecules into chemical pairs.The effectiveness of the model was shown by the formation of two co-crystals of captopril with l-proline and sarcosine.Przybyłek et al. [134] developed a MARSplines algorithm model containing 1D and 2D molecular descriptors and validated it for phenolic acid co-crystals.In a study of Meng-Lund et al. [62], molecular descriptors were computed for six APIs (carvedilol, mebendazole, carbamazepine, furosemide, indomethacin, and simvastatin), as well as for 20 naturally occurring amino acids, to build a partial least squares discriminant analysis (PLS-DA) model.Although this model was capable of accurately predicting 19 out of 20 mebendazole-amino acid co-amorphous combinations, it was limited to the use of amino acids as co-formers.[136] later built a model based on multivariate PLS-DA using the same six drugs and 20 amino acids.The model was then used to predict co-amorphous systems of mebendazole with 29 non-amino acid co-formers where 26 (90%) were correctly predicted.Deng et al. [67] calculated Pearson correlation coefficients and Ridge regressions between several molecular descriptors and co-amorphous forms.It was found that in norfloxacin co-amorphous salt systems, the formation of co-amorphous forms was dependent on the selection of a co-former with more π-conjugated rings, better miscibility with a smaller δ Δ value (ideally less than or close to 5 MPa 0.5 ), and a larger ΔpKa.Wang et al. [61] found that the aminopyridine synthon cannot completely guarantee the cocrystal formation between nimesulide and a series of pyridine analogues.By performing quantitative analysis of Ridge and Lasso regressions, it was discovered that the formation of co-crystals was influenced by several molecular descriptors of the co-formers, ranked in order of their impact as follows: MEPs, h_ema (sum of hydrogen bond acceptor strengths), Kier flex (molecular flexibility), and the horizontal distance of two N atom projections.

Other Approaches
ter Horst et al. [137] proposed a thermodynamic principle-based experimental method to discover co-crystals according to the solubility of their pure components (i.e., starting materials).The method considers the solubilities of individual components to Chambers et al. [136] later built a model based on multivariate PLS-DA using the same six drugs and 20 amino acids.The model was then used to predict co-amorphous systems of mebendazole with 29 non-amino acid co-formers where 26 (90%) were correctly predicted.Deng et al. [67] calculated Pearson correlation coefficients and Ridge regressions between several molecular descriptors and co-amorphous forms.It was found that in norfloxacin co-amorphous salt systems, the formation of co-amorphous forms was dependent on the selection of a co-former with more π-conjugated rings, better miscibility with a smaller ∆δ value (ideally less than or close to 5 MPa 0.5 ), and a larger ∆pKa.Wang et al. [61] found that the aminopyridine synthon cannot completely guarantee the co-crystal formation between nimesulide and a series of pyridine analogues.By performing quantitative analysis of Ridge and Lasso regressions, it was discovered that the formation of co-crystals was influenced by several molecular descriptors of the co-formers, ranked in order of their impact as follows: MEPs, h_ema (sum of hydrogen bond acceptor strengths), Kier flex (molecular flexibility), and the horizontal distance of two N atom projections.

Other Approaches
ter Horst et al. [137] proposed a thermodynamic principle-based experimental method to discover co-crystals according to the solubility of their pure components (i.e., starting materials).The method considers the solubilities of individual components to identify suitable concentration ranges for exploring and discovering potential co-crystals, rather than focusing solely on the specific ratios or proportions of the co-crystal components.The screening steps are shown in Figure 16.Saturation temperature and reference temperature are terms determined from the phase diagram of a co-crystal system at different temperatures in solution.Co-crystals are likely to be formed if the saturation temperature is higher than the reference temperature by >10 • C.This method successfully led to the discovery of new co-crystals, including those involving carbamazepine with isonicotinamide, benzamide, and 3-nitrobenzamide, as well as cinnamic acid with 3-nitrobenzamide.
Pharmaceutics 2023, 15, x FOR PEER REVIEW 27 of 34 reference temperature are terms determined from the phase diagram of a co-crystal system at different temperatures in solution.Co-crystals are likely to be formed if the saturation temperature is higher than the reference temperature by >10 °C.This method successfully led to the discovery of new co-crystals, including those involving carbamazepine with isonicotinamide, benzamide, and 3-nitrobenzamide, as well as cinnamic acid with 3-nitrobenzamide.
In recent years, computational methods have also been developed to screen coformers more accurately and efficiently.A modified, non-bonded interaction energy model was introduced by Deng and co-workers [23] to predict the formation of coamorphous systems.The formation of co-amorphous systems was believed to be driven by the interactions between the API and co-former.Their hypothesis suggests that the strength of these interactions directly influences the likelihood of co-amorphous formation.To quantify these molecular interactions, researchers calculate and denote the non-bonded energy difference (ΔEnon-bonded) between the API and co-former.A higher absolute value of ΔEnon-bonded was indicative of a greater probability of co-amorphous formation.This model was developed based on 105 solvent-based cases involving 13 different drugs, and successfully predicted five new co-amorphous combinations of gefitinib and four new co-amorphous combinations of erlotinib.

Conclusions
Co-crystals, salts and co-amorphous systems are multicomponent solid forms containing one or more types of molecules.They are able to offer significant benefits to the bioavailability, solubility, dissolution, stability, permeability, and other physicochemical and pharmacokinetic properties of APIs.Choosing appropriate coformers for a specific API is essential to develop solid formulations with ideal properties, which generally requires a significant amount of time and resources.Over the past few years, several knowledge-based methodologies, as well as computer-assisted prediction methods, have emerged exhibiting efficacy and accuracy in co-former selection and prediction of the formation of multicomponent solid forms.More strategies have been developed for multicomponent crystalline forms than amorphous forms.Some methodologies, such as HSP, and other new approaches, such as artificial intelligence, are applicable for both solid-state forms.Though none of the single approaches can completely and accurately predict the formation of multicomponent solid forms, each In recent years, computational methods have also been developed to screen co-formers more accurately and efficiently.A modified, non-bonded interaction energy model was introduced by Deng and co-workers [23] to predict the formation of co-amorphous systems.The formation of co-amorphous systems was believed to be driven by the interactions between the API and co-former.Their hypothesis suggests that the strength of these interactions directly influences the likelihood of co-amorphous formation.To quantify these molecular interactions, researchers calculate and denote the non-bonded energy difference (∆E non-bonded ) between the API and co-former.A higher absolute value of ∆E non-bonded was indicative of a greater probability of co-amorphous formation.This model was developed based on 105 solvent-based cases involving 13 different drugs, and successfully predicted five new co-amorphous combinations of gefitinib and four new co-amorphous combinations of erlotinib.

Conclusions
Co-crystals, salts and co-amorphous systems are multicomponent solid forms containing one or more types of molecules.They are able to offer significant benefits to the bioavailability, solubility, dissolution, stability, permeability, and other physicochemical and pharmacokinetic properties of APIs.Choosing appropriate co-formers for a specific API is essential to develop solid formulations with ideal properties, which generally requires a significant amount of time and resources.Over the past few years, several knowledgebased methodologies, as well as computer-assisted prediction methods, have emerged exhibiting efficacy and accuracy in co-former selection and prediction of the formation of multicomponent solid forms.More strategies have been developed for multicomponent crystalline forms than amorphous forms.Some methodologies, such as HSP, and other new approaches, such as artificial intelligence, are applicable for both solid-state forms.Though none of the single approaches can completely and accurately predict the formation of multicomponent solid forms, each approach can offer valuable guidance in selecting the most suitable co-formers.The implementation of multiple screening/prediction approaches can result in considerable enhancements in the effectiveness and accuracy of predicting the formation of co-crystals, salts and co-amorphous systems, as confirmed by previous studies.
It can be expected that in the future, more investigations will continue to promote screening/prediction efficiency by using multiple approaches as well as developing novel models.One direction could be the integration and combination of various approaches mentioned in this review.For instance, combining machine learning algorithms with molecular docking methods could hold potential to improve the precision of multicomponent solid-state forms (co-crystal and co-amorphous) formations.The combination of quantum mechanical simulations with virtual screening techniques offers opportunities to delve deeper into the intricate interplay of thermodynamic stability and intermolecular interactions of multicomponent systems.High-throughput virtual screening methods taking into account advancements in computing power coupled with experimental screening to validate the predictions will speed up the identification of potential co-crystals and co-amorphous forms.
Collaborations between computational scientists and experimentalists are required to arrive at synergistic approaches to co-crystal and co-amorphous discovery.As researchers explore novel combinations of approaches, significant advancements in the field are foreseeable.These collective endeavors are poised to pave the way for the expedited discovery and development of innovative drug formulations, which could have a transformative impact on the pharmaceutical industry.

Figure 1 .
Figure 1.Representation of the different API solid-state forms.

Figure 2 .
Figure 2. Timeline of the development of co-former screening and formation prediction of multicomponent solid forms.

Figure 2 .
Figure 2. Timeline of the development of co-former screening and formation prediction of multicomponent solid forms.

Figure 3 .
Figure 3. Relationship between the relative occurrence of co-crystals (AB neutral (grey)) and salts (A − B + ionic (orange)) and the calculated ∆pK a .Reprinted with permission from ref. [71].Copyright 2012 Royal Society of Chemistry.
. It was found that supramolecular heterosynthons (COOH•••N arom and OH•••N arom ) are preferred over the corresponding supramolecular homosynthons (COOH•••COOH and OH•••OH).Supramolecular synthons enable the prediction of the possible interactions between different molecules and the assessment of their propensity to form stable co-crystals.Through the analysis of complementary synthons in the molecular components, it becomes possible to narrow down the potential co-formers for a given target molecule.Pharmaceutics 2023, 15, x FOR PEER REVIEW 13 of 34 [70] found that robust O-H•••N or N-H•••O hydrogen bonds were observed between pyrimidine N-oxide and carboxylic acids.Therefore, 17 co-former candidates with carboxylic acid groups were chosen and eight aromatic carboxylic acids successfully formed multicomponent crystal forms with minoxidil.

Figure 6 .
Figure 6.MEPs representation displaying charge distribution (red for charge concentration and blue for charge depletion).Reprinted with permission from ref. [87].Copyright 2022 Royal Society of Chemistry.

Figure 6 .
Figure 6.MEPs representation displaying charge distribution (red for charge concentration and blue for charge depletion).Reprinted with permission from ref. [87].Copyright 2022 Royal Society of Chemistry.

Figure 7 .
Figure 7. HBP approach for screening of multicomponent crystal forms.

Pharmaceutics 2023 , 34 Figure 9 .
Figure 9. Co-formers of indomethacin (top) and paracetamol (bottom) co-crystals.Successful and failed cases are highlighted in green and red, respectively.Reprinted with permission from ref. [102].Copyright 2020 American Chemical Society.

Figure 9 .
Figure 9. Co-formers of indomethacin (top) and paracetamol (bottom) co-crystals.Successful and failed cases are highlighted in green and red, respectively.Reprinted with permission from ref. [102].Copyright 2020 American Chemical Society.

Figure 10 .
Figure 10.Three unequal dimensions of a model box.

Figure 10 .
Figure 10.Three unequal dimensions of a model box.

Figure 10 .
Figure 10.Three unequal dimensions of a model box.

Figure 11 .
Figure 11.Multicomponent crystal forms screening of 2,4-dichlorophenoxyacetic acid based on COSMO-RS and MC methods.Reprinted with permission from ref. [106].Copyright 2022 Royal Society of Chemistry.

Figure 11 .
Figure 11.Multicomponent crystal forms screening of 2,4-dichlorophenoxyacetic acid based on COSMO-RS and MC methods.Reprinted with permission from ref. [106].Copyright 2022 Royal Society of Chemistry.

Figure 16 .
Figure 16.Steps for screening the formation of co-crystals based on the solubility of the pure components.Reprinted with permission from ref. [137].*: at constant temperature.Copyright 2009 American Chemical Society.

Figure 16 .
Figure 16.Steps for screening the formation of co-crystals based on the solubility of the pure components.Reprinted with permission from ref. [137].*: at constant temperature.Copyright 2009 American Chemical Society.

Table 1 .
Hydrogen bond based methods for co-former screening and formation prediction of multicomponent solid forms.
•Only 2 out of 75 co-formers (3%) led to the formation of a co-crystal.•Lowsuccessrate was due to artemisinin lacking strong hydrogen bond donors or acceptors.•Structureandintermolecularinteractions of leflunomide were analyzed using Isostar and Mercury to identify favorable functional groups for co-crystallization, resulting in 5 new co-crystals.•CSDmotif search offered 39 potential candidates and MC was used to screen the compatibility.• 13 new co-crystals, 1 salt and 4 co-amorphous systems were identified experimentally.
• FlexCryst software (www.flexcryst.com(accessed on 23 July 2023)) suggests that the prediction of co-crystal formation must satisfy ∆G ≥ −3 kJ/mol for feasibility, instead of ∆G ≥ 0 kJ/mol.•Anextremegradientboostingmodel was developed using 1000 co-crystallization cases and 2083 chemical descriptors.•3newco-crystals of diclofenac were discovered.•A random forest-based co-crystal prediction model was created using a dataset of positive samples from CSD and negative samples from randomly paired molecules.• 2 captopril co-crystals verified the effectiveness of the model.
•Co-crystal formation between nimesulide and pyridine analogues depends on various molecular descriptors of co-formers, with MEP having the greatest impact, followed by h_ema (sum of hydrogen bond acceptor strengths), Kier flex (molecular flexibility), and horizontal distance between two N atom projections.