Kinetics of Lumefantrine Thermal Decomposition Employing Isoconversional Models and Artificial Neural Network

Thermal analysis can be used to determine shelf-life and kinetic parameters in pharmaceutical systems. This work investigates the kinetic of lumefantrine thermal decomposition, an antimalarial, using non-isothermal and isothermal experimental data. The non-isothermal conditions are analyzed applying Vyazovkin method, while isothermal conditions employ models fitting procedure and artificial neural network. Lumefantrine was characterized by powder X-ray diffraction and Fourier transform infrared spectroscopy. The initial stage of lumefantrine thermal decomposition, about 5% of conversion, corresponds to the loss of chlorine and hydroxyl, being correctly predicted by the neural network as a complex event. At room temperature, the D3 model is appropriate to describe the process, once the half-life time is 19 months, in agreement with manufacturer. Isoconversional model determined the activation energy along the whole process while isothermal methodology determined the global value considering the entire process. The results provide important information for the pharmaceutical industry to assay levels of acceptable lumefantrine contents.


Introduction
The investigation of physical and chemical properties as well as the thermal stability of active pharmaceutical ingredients (API) is essential during drug product development. 1,2 Thermogravimetry (TG) and differential scanning calorimetry (DSC) are traditionally used to study drug products, API and excipients concerning their shelflife and kinetic parameters. [3][4][5][6][7][8][9][10] The present work deals with a forced decomposition of the API. To better understand the behavior in shelf time, a long term stabilization study is necessary. A way to assess the long term study is to force the aging of the API. In this analysis, constant heating is assumed to simulate the stability. Naturally, in practice, the API is not submitted to a higher temperature, in storage nor use. However, the kinetic study at high temperatures is used to simulate the necessary aging to have a complete and reliable picture of thermal stability. The stability studies are based on the pharmaceutical system itself, if only the API is studied, the stability is not associated to the formulation. For shelf life, the pharmaceutical system is already present. All kinetic studies deal with higher heating temperatures searching for the pre-determined acceptable minimum amount of API based on the regulatory policies.
Decomposition of solids relates mass loss and time or temperature, occurring in the reactant-product interface with the formation and growth of nuclei. 11 In isothermal decomposition analysis, the kinetic process can be investigated by imperfections on the material surface, where the reactions and nuclei are still preferably formed. Kinetic models, well established in the literature, [12][13][14] are used to explain experimental curves. The usual procedure is to select only one kinetic model to describe the whole event. [12][13][14] In many situations, the use of a single kinetic model can result in unacceptably high residual error. This fact suggests that the single model represents only an approximation of the real phenomenon. 15,16 The models relate fraction of reacted material, normalized to the total mass lost and time, α(t). When the experimental data adjustment of α(t) by only one model is not appropriate, a combination of mechanisms appears as a solution to investigate. For this, a neural network approach is proposed. In this adaptation, the first layer consists of an input neuron receiving experimental data of time for each isotherm and the output neuron determines the degree of conversion. This neural network is adapted to appropriately describe the solid-thermal decomposition with a specific algorithm, which differs from traditional neural network approaches. Our group developed this multilayer perceptron network (MLP) as homemade software in MATLAB ® language and only the interconnection values, called weights, between the intermediate and output layers are updated, which represents the contribution of each kinetic model considered in the network architecture to describe the experimental data of conversion degree in function of time.
The non-isothermal treatment is a reliable method to evaluate the activation energy of thermal decomposition according to the conversion degree. The most commonly used methods are the differential Friedman (FR) method 17 and integral Ozawa-Flynn-Wall (OFW), [18][19][20] Kissinger-Akahira-Sunose (KAS) and Vyazovkin 21 methods. According to International Confederation for Thermal Analysis and Calorimetry (ICTAC) 21 the differential kinetic equation can be treated in the integral form, only if the activation energy and pre-exponential factor are constant during the process. Thus, OFW and KAS are not appropriate for complex processes. The Vyazovkin method is applied at short intervals, allowing an indirect analysis of multiple steps, as the activation energy varies according to α degree. Isoconversional methods require a set of experimental data performed at different heating rates. From this, for some specific conversion, the temperature of the event at these different curves is analyzed. 10 The isoconversinal approach has been used in research areas as environmental chemistry, 22 pharmaceutical industry 23,24 and materials chemistry. 16 Data for degradation processes need appropriate measurement according to ICTAC 21 and kinetic's committee recommendations. [25][26][27][28][29][30][31][32] These analyses allow the evaluation of the material stability, shelf life and consequently the behavior in pharmaceutical formulations.
In this context, we investigated the kinetics of lumefantrine thermal decomposition, employing TG data. Lumefantrine is the first choice for the treatment of uncomplicated and severe Plasmodium falciparum malaria. 33,34 The knowledge of thermal stability is fundamental to properly design the industrial process, since most unit operations require or generate heat, as can occur to obtain tablets containing lumefantrine. 1,2 It should be emphasized there are no data in the scientific literature describing the kinetic of lumefantrine thermal decomposition.

Characterization of lumefantrine API
Powder X-ray diffraction (PXRD) data were obtained using a Shimadzu diffractometer, XRD-7000, Japan, under 40 kV, 30 mA, using Cu Kα, measured from 5-50° 2θ with a step size of 0.02° and a time constant of 1.2 s step -1 , using a graphite monochromator, in a parallel focusing geometry under 30 rpm preventing any remained preferred orientation. Fourier transform infrared spectroscopy (FTIR) analysis was performed in a PerkinElmer spectrometer, Spectrum One, USA, in the range of 600-4000 cm -1 at room temperature with a resolution of 4 cm -1 . The Spartan software 35 was used for IR interpretation and assignments.

Thermal analysis
The DSC curve was obtained in the DSC60 Shimadzu, Japan, under a dynamic nitrogen atmosphere, with a flow rate of 50 mL min -1 , under a heating rate of 10 °C min -1 from 30 up to 400 °C in a closed aluminum crucible. The sample of lumefantrine was about 1.5 mg, accurately weighed. The phenomenon was characterized by the T onset , corresponding to the extrapolated temperature, calculated based on the peak asymmetry of the thermal phenomenon and the ΔH fus corresponds to the involved enthalpy change for the fusion. 36 TG curves were obtained using a Shimadzu DTG60 thermobalance, Japan. Heating rates of 5, 10, 15 and 20 °C min -1 in the temperature range of 30 up to 600 °C were applied for non-isothermal experiments. The specific temperature range for isoconversional treatment was chosen for each curve based in the interval in which only the decomposition process occurs, i.e., the temperature range in which the mass loss was significant.
Isothermal experiments were performed at 255, 260, 267 and 270 °C because at these temperatures one observes a maximum of 2% in the lumefantrine decomposition. This maximum conversion is duly adopted by The International Phamacopeia 37 to study the maximum tolerance of the composition variance in this drug, together with a standard protocol of using heating rate of 10 °C min -1 until reaching the target temperature. Measured isothermal and nonisothermal experiments were under a dynamic nitrogen atmosphere with a flow rate of 50 mL min -1 in an alumina crucible with sample mass range of 2-3 mg, weighted directly in the thermobalance.

Non-isothermal treatment
Solid-state decomposition can be studied assuming the process occurring in first-order kinetic, described by equation 1: where k(T) is the rate constant following the Arrhenius equation and f(α) the reaction model, being α the extent of conversion, defined as, where m 0 = initial sample mass; m f = final sample mass; m t = sample mass at time, t.
The present study employed the integral isoconversional Vyazovkin method. The Vyazovkin nonlinear procedure 11 considers the integral: (3) in which E α is the activation energy at the conversion degree α, R is the gas constant and T the temperature. Together with equation 1 as, (4) in which A is the frequency factor, and assuming the reaction model is independent of the heating rate, β i (i = 1,…,n), we have for a given conversion: Considering A α as a β independent in equation 5, this equation can be represented as a condition of minimum value as (6) The integral of equation 4, the Arrhenius integral, can be expressed as (7) In this equation, x = E/RT and p(x) was calculated using the Senum-Yang approximation 38 of third degree as (8) According to Pérez-Maqueda and Criado, 38 the percentage error of the p(x) function for the 3 rd rational approximation is about 10 -5 percent for these x.

Isothermal treatment
The MLP applied in this study has three layers: one input, one intermediate and one output layer according to the algorithm proposed initially by Sebastião et al. 12,13 The input and output layers have only one artificial neuron, and the intermediate layer has a variable number of artificial neurons, depending on the number of kinetic models to be considered in the process of solids decomposition, in this study, a total of nine models, described in Table 1.
The states of neurons o k , are defined as, with x k being the state of the artificial neuron in the previous layer and w kj the weights between the neurons k and j. The propagation of information occurs, only if an activation function, f, is applied in the artificial neurons. In MLP networks there are also weights called bias, w i0 defined as unity parameters to amplify or reduce the linear correlation of the network. 12,13 For the proposed MLP architecture, the weights w i1 and the bias w i0 were predetermined from the fit of experimental data by each kinetic model in each studied temperature. These weights are the rate constants and linear coefficient, respectively, of the kinetic models. From this, the weights among the input and intermediate layers, do not change in the learning process of the network and can be represented as a weight matrix w 1 as, Also, one has to assume x as the experimental time data for each studied temperature in the following form: (11) From this representation, the matrix product w 1 x determines the state of the artificial neurons in the intermediate layer as (12) These neurons in the intermediate layer should be activated to transfer information to the output layer, analog to the nervous impulse in the brain. Using the kinetic models, represented in Table 1, as the activation functions in the intermediate layer, the neurons in this layer assumes the states, For this proposed methodology, this equation represents the particular case of equation 9 (the general state of a neuron in a network). In equation 13, the terms in the vector (f F1 f F2 f F3 ... f D4 f n ) represent the kinetic models related in Table 1.
The chosen activation functions in neural networks obey some requirements. It is necessary that activation function assumes a predetermined value, generally, f(x) = 0 before calculations of the state of neurons; activate the neurons assuming values near to the unity and be a crescent function , to warranty the energy function optimization. 12 If we choose a linear activation function in the output layer, the network response, given by the state of the neuron in the output layer is (14) being w 2 the interconnection weights vector of the output layer.
The MLP learning process consists of an energy function (E) optimization, from which only the w 2 vector has to be determined. 12 Considering α exp as the experimental data, this energy function is (15) Assuming the w 1 matrix has the kinetic rates determined -ln(1 -α) = kt first-order 0 from the fit of experimental data by the kinetic models and the activation function of each neuron in the intermediate layer as the kinetic models, the w 2 contribution of each kinetic model to describe the whole experimental process can be calculated by minimizing the equation 15, as w 2 B = α exp with B = f(w 1 x). Once the B matrix is an ill-conditioned matrix, the pseudo-inverse can be used to determine the w 2 vector as, with α exp as the experimental conversion degree.

Results and Discussion
PXRD and FTIR spectra of lumefantrine ( Figure 1) confirm the identity of the API used in this study. The characteristic of lumefantrine powder was representative of polycrystalline material, by fully reproducing the X-ray diffraction down to weighted profile R-factor (R wp ) of about 3%. The R wp follows directly from the square root of the quantity minimized, scaled by the weighted intensities. Counting longer increases the mathematical precision in a diffraction measurement. Indeed, as the total number of counts collected for a diffraction pattern is increased, expected R (R exp ) decreases. Paradoxically, counting longer will usually increase the difference between R exp and R wp , thus making worse the model obtained by fitting. This is due to the patterns measured with very large numbers, even minor "imperfections" features that cannot be modeled in the peak shape or peak positions can make it impossible to obtain small χ 2 or R wp values. 40 The characteristic bands of lumefantrine were at 3398 cm -1 , corresponding to O−H stretching, with variable intensity; at 2928 cm -1 , corresponding to C−H linear stretching, with medium intensity; and at 770 cm -1 , corresponding to C−Cl range, with high intensity. At left in Figure 1, the PXRD shows Rietveld fitting graphic (FullProf.2k, Version 6.30, Sep2018-ILL JRC). 41 The fitting well predicts all experimental points. The difference plot is near straight line, stating the material purity.
The DSC curve for lumefantrine presented an endothermic event, not followed by a mass loss in the TG curve, which indicates a physical phenomenon related to melting (T onset = 128.1 °C, ΔH fus = 77.39 J g -1 ), as expected. 37 The enthalpy is calculated based on the endothermic event integration by TA 60 data software version 2.20, Shimadzu, Japan. The decomposition process occurred after this fusion process, starting at about 250 °C ( Figure 2).
The non-isothermal experimental data of α(T) obtained with four different heating rates (5, 10, 15 and 20 °C min 1 ) for lumefantrine are presented in Figure 3. The four curves indicate that the whole decomposition event occurred in only one step. Once the fusion had occurred before the decomposition event, this experimental finding is coherent, indicating mass loss as a continuous process along with the thermal decomposition. All the curves were treated by the non-linear method of Vyazovkin. 11  Figure 4 shows the activation energy determined by the non-linear Vyazovkin method. 11 The residual errors are presented in Table 2. The activation energy is a crescent function along with the complete conversion. This experimental result indicates that as the thermal decomposition process occurs, more energy is necessary to maintain the process. Therefore, any catalyst substance that could decompose lumefantrine is formed. Also, one can note that this result is significant to attest the thermal stability of lumefantrine. For example, to start the thermal decomposition process, it is necessary activation energy of about 105 kJ mol -1 , according to Vyazovkin method, 11 and to maintain the process it is necessary to increase the energy to 110 kJ mol -1 to reach 20% of conversion.
Another non-isothermal analysis was also performed considering a general empirical kinetic model adjustment,     with A as the frequency factor, β for heating rate, E a the activation energy, R the gas constant, T the temperature data and m, n and q as parameters to be determined from an optimization process using experimental data. In this procedure, the frequency factor, activation energy and kinetic model can be determined, considering the activation energy constant during along the whole process. The fit procedure was performed with lumefantrine experimental data at the four heating rates and the obtained parameters are A = 1.308 × 10 4 s -1 , E a = 40.2 kJ mol -1 , m = 0.11876, n = 1.2994 and q = 1.1005, with residual error of 2.57 × 10 -3 . As presented by Cai and Liu 42 and Pérez-Maqueda et al. 43 when the retrieved m, n and q parameters do not correspond to the ideal models, as listed in Table 1, it suggests the process follows a complex kinetic and should not be described as a single ideal model, but preferentially as a combination of kinetic models. From the obtained results, we can note these parameters do not correspond to the parameters of ideal models. So, it was intended to use the neural network methodology to determine the combination of ideal kinetic models that accurately describe the lumefantrine experimental data. Experimental isothermal curves at 255, 260, 267 and 270 °C for lumefantrine are demonstrated in Figure 5. These temperatures are due to the correspondence of the decomposition process beginning, verified by the non-isothermal analysis, corresponding to the maximum of 2% in the decomposition, in full agreement with the acceptable limit for the assay of lumefantrine as described in The International Pharmacopeia. 37 Table 3 contains residual errors of adjustment by the kinetic models and MLP. Based on the found values, the isothermal models should be applied to lower values of conversion rates, as required in the pharmaceutical industry to provide assay levels of acceptable API's contents.
It may be noticed that the kinetic models present small residual errors to fit experimental data, being in the same order of magnitude, i.e., the majority of models presented the same exponent to the power of ten, except the R2 model at 260 and 267 °C and D2 at 270 °C, which presented one order smaller. This result suggests that lumefantrine thermal decomposition process have to be understood as a combined event of several mechanisms. For example, the individual models adjustment at 255 °C (first temperature of isothermal experiments) is presented in Figure 6. From this figure, one can note, the models R2 and R3, together with the Avrami-Erofeev, Am2, and Am4, can explain the experimental data with similar accuracy. Although these ideal models presented similar accuracy to describe the data, they present different physical concept and the neural network can measure the contribution of each kinetic model to describe the whole process, considering all of them in the event. Note that the time of exposure to the heating is the determining factor for thermal decomposition in isothermal mode. It took 66 min for degradation of approximately 30%, as an average value, for all four temperatures (255, 260,  267 and 270 °C). In contrast to isoconversional mode, the limiting factor is the heating rate. All these models were considered in intermediate layer of the proposed MLP, with w 1 matrix determined by the rate constant, k, and linear coefficient, k0, of kinetic models by fitting experimental data in each isothermal curve. The activation function of each neuron in this layer corresponds to the kinetic model, applied in the w 1 x term. The w 2 vector was determined as described in equation 16, and this term also corresponds to the kinetic model corrections by the neural network. The kinetic model improvement by the network takes into account the residue of the sample, basically a charcoal product. Once the asymptotic value of the models is unitary, the w 2 vector can correct this value and additionally provide the contribution of each model to describe the experimental data.
The corrected rate constants of the w 1 matrix were used to calculate the activation energy of lumefantrine thermal decomposition in this initial stage, values presented in Table 4. It is important to emphasize this value refers to the activation energy necessary to initiate lumefantrine decomposition. According to isoconversional experimental data, the temperature for 5% of conversion is 274.6 °C. Then, we can consider the process, occurring at the range of 25 up to 255 °C, with activation energy relatively constant for any conversion less than 2% (the limit for the assay of lumefantrine as described for World Health Organization). 37 These values showed in Table 2 are much smaller than the activation energy determined for conversion higher than 5%, justifying the approximation. Also, the activation energy determined from the Cai and Liu 42 method is about 45 kJ mol -1 , which is in agreement with the neural network results. Taking into account those activation energy values, the rate constants and consequently the half-life time, can be determined for the process at 25 and 255 °C, using Arrhenius equation. Thermal stability can also be demonstrated by calculating a prediction of the API half-life time reaction rate at, e.g. 25 °C, based on the determined kinetic parameters. This temperature was chosen to attend the United States Pharmacopeia for API storage, which states the temperature between 15-30 °C. 44 From the determined kinetic parameters at 255 °C, Table 4 presents the half-life time of lumefantrine at both temperatures. At 25 °C, it seems D3 is the more appropriate model to describe the decomposition process, once the half-life time calculated is 19 months, in agreement with the period established by the manufacturer.
Individually analyzing the kinetic models, the R2 and R3 models adjusted the phenomena better for lower temperatures and D2 (diffusion) and R3 models for higher temperatures. From these results, we can note for room temperature, the diffusion mechanism is a rate-limiting process, once the system is in a solid phase. As the fusion temperature is reached, the mechanisms change to R2 and R3 model. After 270 °C, the diffusion mechanism is again the dominant process, probably due to some intense interaction energy to be overcome to continue the decomposition process. The Rn models are mathematical functions to describe reactions in the boundary phase, occurring when the diffusion process is extremely fast, which prevents an adequate interaction between reactants in the interface. For this, the process is controlled by the rearrangement of the reactants in the limitrophe phase. These models are known as geometrical contraction models and the functions are similar to the models described by reaction order mechanisms, applied to a homogeneous system. So, as the lumefantrine system melts at 128 °C, the accuracy of R2 and R3 models to describe the experimental data at higher temperatures is related not to the physical meaning of R2 and R3 models, but to the similarity of these models with the reaction order mechanisms used to describe homogeneous systems. This work aims to evaluate the models as mathematical functions to accurately fit the experimental data, providing accurate rate constants to be applied at storage temperature.
The contribution of the models for lumefantrine isothermal decomposition process for each temperature (normalized values) indicates that at 255, 260 and 267 °C the most significant contribution was the R2 model. At 270 °C, the most contributing mechanisms were the D2, diffusion, and R3 models (Figure 7). At this point, it is essential to observe the initial activation energy was found to be very low, 10-20 kJ mol -1 , but reached about 100 kJ mol -1 at the conversion of 0.05 and then increased with the extent of conversion after that. The behavior seems to indicate some elements of diffusion or mass transport rate-limited degradation, maybe due to the leaving of cleaved fragments which initially occurs at the melted surface and moves to the deeper inside of the sample.
Once the sample is liquid at temperatures higher than 255 °C, the R2 and R3 models have to be considered as reaction order kinetic models for homogeneous systems. When the temperature gets higher, the chemical reaction rate increases, and now the removal of the product becomes a rate-limiting step in the process. Also, larger size sample study was performed and at 255 °C, a mass of 3.570 mg (43% larger) was tested. For this sample, the neural network results pointed out a better fit of R2 model (smaller residual error of 0.154) with a higher contribution of this model in the network to describe the process (30%), indicating the process indeed occurs following a second order homogeneous kinetic instead of the diffusion mechanism.
The observation of the residual errors and contributions of each kinetic model suggests the R2 and R3 models, as reaction order mechanisms, are preferable to describe the process before 270 °C, where the diffusion models became more expressive. For the decomposition pathway evaluation, the sample was thermally treated at 314, 318 and 324 °C in separate experiments. The FTIR spectra of lumefantrine, in each temperature ( Figure S1, Supplementary Information), corroborate with the fitted kinetic model, since in the initial stage of the decomposition, only the loss of chlorine and hydroxyl group was verified. The reduction in the bands at  3394 cm -1 (C−OH) and 800-600 cm -1 range (C−Cl), supports the conclusion. Based on the theoretical FTIR obtained in Spartan simulation software, 35 the OH and the Cl groups are the lost groups at the decomposition starting point. The loss of these groups promotes a contraction in the area and volume of the remaining material. The Molispiration platform 45,46 was also used to determine the volumes of structures before and after the beginning of its thermal decomposition. The initial volume was 455.61 Å 3 , and the volume after chlorine or hydroxyl loss was 434.03 Å 3 , whose volumes were 24.24 and 18.01 Å 3 , respectively. The structures are shown in Figure 8.
It is important to emphasize that our intention is not to generalize the volume contraction behavior with cleavage of any small molecule in the decomposition process. In the present work, the mass loss comprises one Cl (chlorine) and OH (hydroxyl) group, the total amount of about 51-52 g mol -1 in a total molecular mass of about 528.94 g mol -1 . It is about 10% of the mass loss, verified in the TG experimental curves. For this case, the volume shrinkage due to the cleavage seems reasonable, as the phenomenon is also detectable in simulated IR using Spartan 14 simulation software. 35 The theoretical IR is compatible with the loss of Cl and OH, resulting in a volume change, at least for this API. The R2 and R3 models are for reactions controlled by phase boundary, which is more typical for solid-state reactions. In this case, they presented a better adjustment of experimental data in smaller temperatures, being considered as reaction order mechanisms, both for the neural network and individually. Although the system is not in the solid phase at this temperature, the models agree with the physical phenomenon of observable Cl and OH loss, associated with final total volume reduction, for the studied system.

Conclusions
The kinetics of lumefantrine thermal decomposition employing isoconversional model and the artificial neural network was defined. The initial stage of the thermal decomposition pathway was studied using isothermal experiments, which suggests a loss of chlorine-3 and hydroxyl group. A general non-isothermal empirical method and neural network results suggest the process occurs as a complex event. From the neural network results, the R2 and R3 models occur at smaller temperatures and diffusion and R3 models at higher temperatures, as a combined event. These models are considered as preferable ones to describe the process at the studied temperatures. It was possible to know the activation energy to initiate the process of lumefantrine thermal degradation by the isothermal conditions as well as the activation energy at 20-80% of conversion by the isoconversional method. The isothermal method is suitable to provide information about the lumefantrine acceptable limit assay (lower conversion rates), being the isoconversional methods applicable for higher conversion rates.

Supplementary Information
The FTIR spectra of lumefantrine, in each temperature, is available free of charge at http://jbcs.sbq.org.br as PDF file.