Process Optimization, Genetic Programming Modeling and Reliability Assessment of Bioactive Extracts Recovery From Phyllantus Emblica

1.0. Introduction

Phyllanthus emblica, also known as Emblica officinalis or Indian gooseberry, is a plant widely recognized for its medicinal properties. P. emblica is native to tropical south-eastern Asia and it is one of the most extensively studied plants. All parts of P. emblica, including its fruits, have medicinal effects and are used to treat diarrhea, jaundice, inflammation, and act as a strong rasayana (techniques for lengthening lifespans) in Ayurveda [1]. There are documented reports of the applications of P. emblica in memory enhancement, respiratory, skin and ophthalmic disorders, and detoxification, including snake venom [2]. P. emblica holds a unique position in Indian (Ayurvedic), Turkish, Unani, and Tibetan medicinal systems and has therapeutic and healing potentials, making it a valid research option for the development of novel drug formulations with low side effects [1]. The plant is widely used for its nutritional and therapeutic properties, including antioxidant, antidiabetic, and antimicrobial effects [3]. However, the diverse activities of the polyphenols that are present in medicinal plants to a greater extent are responsible for their disease ameliorating properties [4] and these make them essential components of bioactive extracts. Phenolic compounds of plant origins have gained significant attention due to their numerous health-promoting properties and potential therapeutic benefits [5,6]. Hence there arises the need for the recovery of high quality of this class of bioactive compounds from P. emblica which could find applications in diverse fields of endeavor.

Heat-assisted extraction (HAE) technology has been a subject of interest and is among the commonly used conventional technologies in the extraction of bioactive compounds from various plant sources [7,8,9,10]. HAE technology (HAE-T) utilizes heat to enhance the extraction efficiency and improve the yield of desired compounds. Other frequently employed traditional techniques include maceration, Soxhlet extraction, and hydro-distillation [8]. Several parameters significantly influence the HAE efficiency in terms of yield and quality. These parameters include temperature, time, particle size, solid to liquid ratio, solvent type and extraction cycle [11]. Polar solvents such as water, methanol, ethanol, and acetone are commonly used for extracting polar bioactive compounds like phenolics and flavonoids [12]. On the other hand, non-polar solvents like hexane, ethyl acetate, and chloroform are preferred for extracting non-polar compounds such as essential oils and lipids [13]. In many cases, a single solvent may not be sufficient to extract the full range of bioactive compounds from a complex matrix. Therefore, solvent mixtures are often employed to improve extraction efficiency and broaden the spectrum of extracted compounds. Commonly used mixtures include ethanol-water, methanol-water, and ethyl acetate-methanol-water. These mixtures take advantage of the varying solvent polarities to extract a wider range of bioactive compounds with different polarities.

Literature searches have revealed several documented works on the use of HAE for the recovery of bioactive compounds from different plant sources. For instance, phenolic-rich extracts have been successfully recovered from Carica papaya leaves, Azadirachta indica leaves, Parkia speciosa pod Senna alata and Enantia chlorantha stem bark [4,9,14,15] using the HAE-T. HAE-T offer certain advantages, such as improved extraction efficiency and higher yields, however, may require longer extraction times, higher energy demands, and increased solvent requirements [14]. The HAE-T is notwithstanding, of interest from industrial-point of view because of its relative scalability, vessel availability and ease of maintenance [16]. HAE-T is among the most affordable extraction technologies in developing countries because of its low initial purchase cost and simplicity of operation. Hence there is a need to optimize the HAE process, enhance the yield of bioactive compounds, and improve the overall efficiency of bioactive extract recovery from P. emblica and other plant sources. Proper control and optimization of the HAE parameters are necessary to maximize the extraction efficiency and overcome the potential draw backs associated with this technology.

Response Surface Methodology (RSM) has been widely applied in various fields to optimize processes and reduce experimentation costs. RSM, coupled with factorial designs, offers a statistical approach to efficiently explore and optimize multiple process variables [17]. RSM is an effective approach for optimizing processes when multiple factors and interactions influence the desired outcomes [18,19]. One of the primary advantages of RSM is its ability to generate valuable insights from a relatively small number of experiments, providing significant information with reduced costs compared to traditional one-variable-at-a-time optimization methods [20]. RSM has found applications in areas such as solvent extraction, engineering optimization, machining process parameters, fermentation processes, and hydrogen production from ethanol thermochemical conversion [21,22]. RSM also aids in identifying optimal conditions and understanding the relationships between process variables and responses, facilitating process optimization and modelling of complex systems [23]. RSM experimental designs such as Box-Behnken design have ability to efficiently identify the curvature of the response surface and provide robustness against potential errors in the estimation of the model parameters. Its balance between accuracy and cost-effectiveness makes it a popular choice in designing experiments and optimizing process parameters [9].

Process modelling helps engineers design and analyze complex systems, improve efficiency, reduce costs, optimize resource allocation, and make informed decisions during the design and operation of engineering processes. Multi-gene genetic programming (MGGP) is one of the promising machine learning techniques that are available for creating mathematical or computational representations of physical processes or systems. MGGP is a variant of genetic programming that offers promising capabilities for behavioral modeling and engineering system analysis. MGGP combines the model structure selection ability of standard genetic programming with the parameter estimation power of classical regression to capture nonlinear relationships between input and output variables [24,25]. It effectively searches for suitable relationships between input and output datasets, independent of the physical background of the data [26] and generates mathematical models of predictor response data that are "multi-gene" in nature, involving linear combinations of low-order non-linear transformations of input variables[25]. MGGP has been successfully applied to various engineering problems, such as behavioral modeling of structural engineering systems [24], identifying the dynamic prediction model for an overhead crane [27], and short-term load forecasting in power systems [28]. However, works on the application of MGGP in modelling extraction processes are very scarce in the literature.

Predictive models are widely used in various fields to make informed decisions and forecasts based on available data. However, assessing the reliability of a predictive model is crucial to ensure its accuracy, robustness, and generalizability. Reliability assessment involves evaluating the performance, stability, and limitations of the model [29]. One of the popular techniques of assessing predictive model reliability and robustness is via Monte Carlo simulation. Monte Carlo simulation involves randomly sampling values from probability distributions of input variables to generate multiple simulated scenarios. This technique provides a comprehensive assessment of the reliability and uncertainty of predictive models by considering the variability in input variables. It also allows for a thorough understanding of the model’s behavior under different scenarios and helps in making informed decisions based on the range of predicted outcomes. Monte Carlo simulation is widely used in various fields, including finance, engineering, physics, economics, and many others. Its ability to handle complexity and uncertainty makes it an invaluable tool for assessing the reliability and robustness of models, allowing decision-makers to make more informed choices and manage risks effectively. This technique has been found useful in the assessment of the robustness and sensitivity evaluation of polymer composite BBD-RSM predictive models [4], GPR and ANN compressive strength predictive models [30], biodiesel yield CCD-RSM predictive model [16] and techno-economic MGGP predictive models [30].

Despite interesting bioactivities displayed by P. emblica leaf extracts in the management of broad scope of ailments, the protocol for its bioactive recovery which is necessary for product standardization has not been reported. Also, the engineering works which include process parameter optimization, process modelling and prediction; and model reliability assessment of the extract recovery process are still not available in the literature at present. These engineering works, however, form the basis for the design, optimization, scale-up, and techno-economic feasibility analysis of an integrated plant for commercial production of P. emblica extracts and which are presently lacking to the best of authors’ knowledge. Therefore, the objectives of this work are to (i) investigate the effects of process parameters on the recovery of bioactive extracts from P. emblica leaves using the Box-Benhken design (ii) optimize the bioactive recovery process with the assistance of desirability function algorithm and characterize the bioactive extract with HPLC (iii) model the recovery process with BBD-RSM and MGGP (iv) conduct a comprehensive model reliability assessment and process variable sensitivity analysis of the better predictive model between BBD-RSM and MGGP.

2. Materials and Methods

2.1. Materials and Chemical Reagents

Fresh leaves of P. emblica were harvested from a medicinal farm plantation and were rinsed with running water in order to remove all the dirt when received in the laboratory. The P. emblica leaves were afterward spread on a tray in the laboratory for about 2 weeks to allow them dry by natural convection. The dried P. emblica leaves were latter ground and sieves to particle size of 0.1 mm and the screened samples were kept inside a black polythene bag before experimentation. All the chemicals used in the laboratory, including absolute ethanol (99%), sodium carbonate, (Na₂CO₃), ascorbic acid, Folin Ciocateu phenol reagent, glacial acetic acid, 2, 4, 6-tripyridyl-s-triazine, sodium acetate trihydrate, HCl, and FeCl₃.6H₂O, were of analytical grade and were purchased from Sigma Aldrich located in Poole, England. Additionally, the HPLC standards employed in the study consisted of betulinic acid, gallic acid, caffeic acid, ellagic acid, chlorogenic acid, ferulic acid, rutin, and quercetin. Distilled water was used though out the laboratory investigations.

2.2. Preliminary Selection of Extraction Solvent

The selection of the best solvent for the heat-assisted extraction (HAE) of antioxidant phytochemicals from P. emblica leaf sample was investigated here. The two extreme pure solvents that were considered were the absolute ethanol (99%) and 100% distilled water. Other solvent mixtures investigated included the 20%, 40%, 50% and 80% ethanol solutions. The heat-assisted extraction of bioactive compounds from P. emblica leaves was conducted at fixed HAE parameters of 40 ^oC operating temperature (OT), 45 min extraction time (ET) and solid to liquid ratio (S:L) of 1:20 g/mL. P. emblica leaf particle size of 0.1 mm was used for all investigations. The responses selected for these investigations were the bioactive extract yield (EY), total phenolic content (TPC) and antioxidant activity (AA) of the extracts. The procedures for the HAE and EY; and TPC and AA quantifications are detailed in Section 2.3 and Sections 2.6 - 2.8 respectively. All investigations were carried out in duplicates and the results were reported as mean ± standard deviation.

2.3. Heat-assisted Extraction of Bioactive Compounds

In the laboratory, the heat-assisted extraction (HAE) of bioactive compounds from P. emblica leaves followed a procedure similar to the method described by Adeyi et al. [16], utilizing a water bath. The process involved placing a known quantity of P. emblica leaf sample and solvent into a beaker and subjecting it to elevated temperature for a specific duration while continuously stirring. To prevent solvent loss, the beaker was covered with foil during the experiment. The HAE experiments were conducted at OT of 40 ^oC, ET of 45 min and S:L of 1:20 g/mL as described in Section 2.2 for the preliminary selection of the best HAE solvent/mixture for the recovery of antioxidant phytochemicals from P. emblica leaves. However, for the study of the interactive effects of the HAE parameters on the extractability of bioactive antioxidants from P. emblica leaves and subsequent process optimization, the solvent used was the best mixture selected according to Section 2.2 and each extraction experiment was conducted based on the experimental conditions outlined in the Box-Benhken design, as explained in Section 2.4. After the designated residence time, the P. emblica leaves fiber was separated from the extracts through centrifugation (at 25 ^oC, 500 rpm, for 10 min). The resulting P. emblica leaves extract was then poured off and stored in plastic bottles for subsequent analysis.

2.4. Box-Behnken Experimental Design and Modeling

The design and investigation of HAE process for extracting bioactive compounds from P. emblica leaves were conducted using the Box-Behnken design-response surface methodology (BBD-RSM). The HAE process involved three key parameters: operating temperature (OT/X₁), solid-to-liquid ratio (S: L/X₂), and extraction time (ET/X₃). These parameters were studied at three levels (-1, 0, +1) in the BBD-RSM. The specific values for OT were 35, 45, and 55 ^oC, for S: L were 20, 35, and 50 mL, and for ET were 100, 150, and 200 min. These levels served as the independent HAE parameters in the study. The measured responses included the total phenolic content (TPC, mg GAE/g d.w), extract antioxidant activity (AA, μM AAE/g d.w), and extract yield (EY, %). The relationships between the HAE parameters and the responses were defined using Eq. (1), where the response parameter is denoted as Y.

$$\mathrm{Y}={\mathrm{b}}_{\mathrm{o}}+\sum _{\mathrm{i}=1}^3{\mathrm{b}}_{\mathrm{i}}{\mathrm{X}}_{\mathrm{i}}+\sum \sum _{\mathrm{i}<\mathrm{j}=1}^3{\mathrm{b}}_{\mathrm{i}\mathrm{j}}{\mathrm{X}}_{\mathrm{i}}{\mathrm{X}}_{\mathrm{j}}+\sum _{\mathrm{i}=\mathrm{j}}^3{\mathrm{b}}_{\mathrm{i}\mathrm{i}}{\mathrm{X}}_{\mathrm{i}}^2$$

(1)

The coefficients $b_o$, $b_i$, $b_{ij}$ and $b_{ii}$ represent the linear, cross-product, and quadratic model coefficients. The HAE parameters in coded form are represented as $X_i$ and $X_j$.

The statistical significance of all model coefficients was determined at a 95% confidence level using the Design-Expert software. The adequacy of the models was evaluated using ${\mathrm{R}}^2$, ${\mathrm{R}}_{\mathrm{a}\mathrm{d}\mathrm{j}}^2$ and Pred ${\mathrm{R}}^2$ as measures of model performance.

2.5. MGGP-Based Predictive Modeling Theory and Model Development

2.5.1. MGGP Theory

Multi-gene genetic programming modelling (MGGP) is a computational approach that combines multiple genetic programming sub-models to evolve complex models capable of capturing the intricacies of the system under investigation. The theory behind MGGP encompasses various aspects, including the representation of individuals, the evolution operators, and the fitness evaluation [4].

One of the key elements in MGGP is the representation of individuals. In MGGP, each individual consists of a collection of sub-models, also known as genes, which are combined to form the overall model. Each gene represents a different aspect of the system and contributes to the overall functionality and performance of the model. The combination of these genes allows MGGP to capture the interactions and dependencies among various variables in the system. Figure 1 illustrates a typical MGGP gene structure.

The model structure consists of independent variables denoted as X₁, X₂, and X₃, as shown in Figure 1. These variables are linked together using mathematical node functions such as subtraction, cosine, tangent, and addition. The mathematical representation of an MGGP regression model is given by Eq. (2).

Y = X₀ + X₁ × tree1 + … + x_k × treeK

(2)

In the equation, the constant term or bias is denoted as x₀, while x₁,..., x_M represent the gene coefficients or weights. The value of M corresponds to the number of genes present in the current or i^th set of inputs.

The evolution operators in MGGP play a crucial role in the creation and modification of individuals over successive generations. These operators include crossover and mutation. Crossover involves the exchange of genetic material between individuals, specifically the sub-models or genes, to create offspring with a combination of characteristics from the parent individuals. Mutation introduces small random changes in the sub-models or genes, allowing for exploration of the search space and potentially discovering new and improved models. Fitness evaluation is an essential component in MGGP to guide the evolution process. The fitness function assesses the performance of individuals based on their ability to solve the problem at hand. It is typically defined based on the specific objectives of the modelling task, such as accuracy, prediction error, or fitness criteria specific to the problem domain. Individuals with higher fitness values are more likely to be selected for reproduction and contribute to the next generation. The overall goal of MGGP is to iteratively evolve individuals with improved fitness through the application of evolution operators and fitness evaluation. This iterative process aims to search for optimal or near-optimal models that provide accurate representations of the system under study.

2.5.2. MGGP Model Development

A Genetic Programming Toolbox for Identification of Physical Systems (GPTIPs), which is open-source software, was utilized in Matlab R2014b® to construct an appropriate MGGP model that captures the relationships within the observed data. The configuration and settings for the modelling parameters of MGGP employed in this study were optimized for the prediction and are as detailed in Table 1. In the course of model development, the experimental data sets (each for TPC, EY and AA) were divided into training (50%) and testing (50%) data subsets.

The best model was selected on the basis of providing minimum fitness value on the training data and developed model simplicity. The effectiveness of the mathematical models developed using MGGP was assessed by calculating the coefficient of determination (R²) and root mean square error (RMSE) values as indicated in Eq. (3) and Eq. (4) respectively.

$${\mathrm{R}}^2=1-\frac{{\sum }_{i=1}^n\left(xi-y_i\right)2}{{\sum }_{i=1}^n{\left(x_i-x\right)}^2}$$

(3)

$$\mathrm{RMSE}=\sqrt{\frac{1}{n}{\sum }_{i=1}^n{\left(x_i-y_i\right)}^2}$$

(4)

where, $x_i$, $y_i$, x and n are observed, predicted output, mean of the observed and number of observations respectively.

2.6. Numerical Process Optimization and Validation

The optimal HAE conditions for recovering bioactive antioxidants from P. emblica leaves were determined using the numerical optimizer provided by the Design-Expert software. The desirability function algorithm, following the method of Adeyi et al., [4], was utilized for this purpose. The optimization criteria involved manually selecting the optimization goal, weight, and importance for both the dependent and independent variables.

In order to find an optimum solution, the goals are combined into an overall composite function D (x), called desirability function. The value of desirability function for n responses can be defined as:

$$D\left(x\right)=(d_1\times d_2\times d_3\dots \dots \dots \dots \dots \dots \dots \dots \dots .d_n{)}^{1/n}$$

(5)

For any response, the desirability varied from zero to one. A value of one represents the ideal case, while zero indicates that one or more responses fall outside the desirable limits.

In this study, the goal was to simultaneously maximize the total phenolic content (TPC), extract antioxidant activity (AA), and extract yield (EY). To validate the predicted optimal HAE conditions, laboratory experiments were conducted. The predicted results were compared to the validated data using the relative standard deviation (RSD) defined by Eq. (6).

$$\%\mathrm{R}\mathrm{S}\mathrm{D}=\frac{{SD}_{1-2}}{M_v}\ast 100$$

(6)

The RSD calculates the standard deviation (${SD}_{1-2}$) and determines the similarity or dissimilarity between the predicted and validated data by comparing the mean value ($M_v$) obtained from both datasets.

2.7. Determination of Extract Yield

The determination of the extract yield (EY) followed the procedure outlined by Alara et al. [22]. In summary, the decanted extracts obtained from the separation of P. emblica leaves fiber and extracts in Section 2.3 were collected. These extracts were then concentrated and dried in a convective oven dryer until a constant weight of P. emblica leaves extracts was achieved. The EY (%) was calculated using Eq. (7).

$$\%\mathrm{E}\mathrm{Y}=\frac{W_1}{W_2}\ast 100$$

(7)

In Eq. (7), the EY represent the percentage yield of P. emblica leaves extracts, $W_1$ refers to the weight of P. emblica leaves extracts (in grams), and $W_2$ represents the weight of P. emblica leaves (in grams) used in the HAE process.

2.8. Determination Total Phenolic Content

To determine the TPC of aqueous extracts, the Folin-Ciocalteu method was used following [31]. To a tube containing 1 mL of ECBEs, which had been diluted with distilled water at a ratio of 1:10, approximately 1.8 mL of Folin-Ciocalteu reagent was added, and the mixture was allowed to stand for 5 minutes. Then, 1.2 mL of 7.5% (w/v) sodium carbonate was added, and the mixture was stirred and incubated at room temperature of 20 ^oC for 1 hour. The TPC (mg GAE/g d.w) was measured spectrophotometrically at an absorbance of 765 nm and calculated using the pre-constructed gallic acid curve (y = 0.037x + 0.0025, R² = 0.997) according to Eq. (8).

$$\mathrm{T}\mathrm{P}\mathrm{C}=\frac{\mathrm{C}\ast \mathrm{V}}{\mathrm{M}}$$

(8)

In this equation, C (mg/mL) represents the concentration of gallic acid as interpreted from the calibration curve, V (mL) is the volume of the solvent, and m (g) is the weight of the dried P. emblica leaves used in the laboratory experimentation.

2.9. Determination of Antioxidant Activity

To evaluate the antioxidant activity of extracts, the ferric ion reducing antioxidant power (FRAP) method was used, as described by Uddin et al. (2016). The FRAP methodology is based on the ability of an antioxidant to reduce a ferric ion, and was used in this study to measure the reducing capacity of P. emblica leaf extracts. To carry out the FRAP assay, a FRAP reagent was prepared by mixing 0.1 M acetate buffer, 0.01 M TPTZ solution, and 0.02 M FeCl₃*6H₂O solution in a 10:1:1 ratio. Then, 50 µL of extract (diluted with distilled water at 1:10) was mixed with 1.5 mL of the FRAP reagent to create a reaction mixture. The mixture was incubated for 30 minutes at 25 °C, and then the absorbance of the reaction mixture was measured at 593 nm using a UV-Vis spectrophotometer. To determine the AA (expressed as µM AAE/g d.w.), a calibration curve of ascorbic acid standard was used, with equation y = 0.0051x + 0.0015 and R² = 0.995.

2.10. High Performance Liquid Chromatography Profiling of Extract

For the characterization of P. emblica leaves extracts, HPLC analysis was performed following the procedure outlined by Mirunalini [3]. The HPLC analysis was conducted using the bioactive extract obtained under optimized process conditions. Prior to the analysis, the extract was filtered using a 0.45 µm membrane to remove any unwanted particulate matter, such as fiber particles. The filtered extract (approximately 10 µL) was then injected into an HPLC system equipped with a UV diode array detector (UV-DAD) set to a wavelength range of 190-800 nm. The HPLC system, a Shimadzu Corporation instrument from Kyoto, Japan, utilized a VP-ODS column (150 × 4.6 mm, 5 µm particle size) to separate the bioactive compounds. The detection of these compounds occurred at a wavelength of 220 nm. The elution of the compounds from the column was achieved using a mobile phase consisting of 0.2% v/v formic acid and acetonitrile. During the analysis, the resulting chromatogram was recorded. To process the HPLC data, Windows LC solution software was utilized. This software facilitated the identification and quantification of the bioactive compounds present in the extract. The retention times and peak areas of the compounds were compared to those of known reference materials to verify their identity. Furthermore, a chromatographic fingerprint of the bioactive plant extract was generated by plotting the retention times and peak areas of the identified compounds.

2.11. Reliability Assessment of BBD-RSM-based Predictive Models

The BBD-RSM-based models were selected as marginally better models for the prediction of bioactive extract recovery from P. emblica leaves based on the statistical evaluation parameters and hence their reliability and robustness will be investigated in this section. Therefore, the reliability of the predictive equations developed using the BBD-RSM approach to estimate the TPC (mg GAE/g d.w), AA (μM AAE/g d.w) and EY (%) of P. emblica leaves extract, was assessed through Monte Carlo simulations performed in the Oracle Crystal Ball software (OCBS). Additionally, the sensitivity of the models to each process variable was examined. In the OCBS, the uncertainty and sensitivity analyses were conducted by assigning the input variables, including ET (min), OT (^oC), and S/L (g/mL), to the assumption cells. On the other hand, the forecast cells were set as TPC (mg GAE/g d.w), AA (μM AAE/g d.w) and EY (%). Uniform distributions were assumed for all the input variables, and a total of approximately 30,000 iterations were executed to ensure a low mean standard error in TPC (mg GAE/g d.w), AA (μM AAE/g d.w) and EY (%).

3. Results

3.1. Preliminary solvent mixture selection for the bioactive compound’s extraction

The selection of ethanol as a co-solvent for recovering bioactive extracts from P. emblica leaves stems from its versatility and effectiveness in extracting bioactive compounds from medicinal plants. These advantages include cost-effectiveness, safety, efficiency, preservation of bioactivity, scalability, eco-friendliness, and regulatory approval, among others. Table 2 presents the relative ability of different solvent mixtures in the recovery of bioactive compounds from P. emblica leaves. As shown in the table, all ethanol-water mixtures including the 100% absolute ethanol and distilled water were able to recover bioactive antioxidants from P. emblica leaves. This shows that a good number of bioactive compounds in P. emblica leaves are polar in nature. The mean values of the extract yield (EY), Total phenolic content (TPC) and Antioxidant activity (AA) were in the range 5.87 – 12.07 %, 13.67 – 24.25 mg GAE/g d.w and 1.98 – 3.76 μM AAE/g d.w respectively for investigations conducted at OT = 40 ^oC, ET = 45 min and S:L = 1:20 g/mL. Close examination of Table 1 revealed that the bioactive extract obtained with absolute ethanol as solvent possessed the least EY, TPC and AA is (EY = 5.87%, TPC = 13.67 mg GAE/g d.w, AA = 1.98 μM AAE/g d.w) while the bioactive extract recovered with 50% ethanol-water mixture demonstrated highest bioactive EY, TPC and AA (EY = 12.07%, TPC = 24.25 mg GAE/g d.w, AA = 3.76 μM AAE/g d.w). Although, the 100% distilled water (EY = 8.33%, TPC = 19.11 mg GAE/g d.w, AA = 2.81 μM AAE/g d.w) did not perform very badly, it was still far below the 80% ethanol-water (EY = 8.41%, TPC = 20.10 mg GAE/g d.w, AA = 2.99 μM AAE/g d.w), 60% ethanol-water (EY = 11.36%, TPC = 23.76 mg GAE/g d.w, AA = 3.66 μM AAE/g d.w), 40% ethanol-water (EY = 11.11%, TPC = 22.98 mg GAE/g d.w, AA = 3.52 μM AAE/g d.w) and 20% ethanol-water (EY = 8.67%, TPC = 19.73 mg GAE/g d.w, AA = 3.31 μM AAE/g d.w) solvent mixtures in the recovery of bioactive extract from P. emblica leaves.

These results clearly demonstrated varied solubility of bioactive compounds, TPC and AA in different ethanol concentrations (40 - 60% ethanol solutions). This work is in close agreement with previous works [32]. Huaman-Castilla et al., [32] pointed out that solvent polarity and ability of solvents to form hydrogen bonds with plant metabolite significantly impact the solvation capacity of solvents which consequently determines the extractability of polyphenols in plant matrices. This present study has however, shown that ethanol concentrations within 40 - 60% range enhanced the extraction of phenolic compounds due to the optimal solubility and stability of these compounds compared with the solvent extremes (absolute ethanol and 100% distilled water). Although the interesting range of ethanol-water solvent concentration for high bioactive solubility is 40 – 60%, the 50% ethanol-water solvent mixture demonstrated the best phytochemical extractability from P. emblica leaves. However, the literature has documented various optimal ethanol-mixture concentrations for the recovery of bioactive extracts from different medicinal origins. For example, a study conducted by Altıok et al. [33] on the extraction of bioactive compounds from olive oil demonstrated that a 70% ethanol solution optimized the extractability of total polyphenols, resulting in extracts with the highest antioxidant capacity. In the same vein, Cheaib et al. [34] concluded that 50% ethanol solution was best for the recovery of bioactive compounds from apricot pomace. Regardless, the addition of ethanol to water as solvent mixture greatly improved the polyphenol-rich antioxidant extraction from P. emblica leaves. This observation was attributed to the impact of ethanol on cell permeability, brought about by alterations in the phospholipid bilayer of the cell membranes, resulting in both chemical and biophysical modifications to the cell membrane [34]. Therefore, the 50% ethanol-water solvent mixture was hence forth used for the selection of the best set of process parameters for the recovery of bioactive compounds from P. emblica leaves in the subsequent Section.

3.2. Effects of process parameters on the recovery of bioactive extract form P. emblica leaves

The BBD-RSM was used to investigate the effects of process variables on the process responses. The variables selected for this investigation were the operating temperature (OT, solid to liquid ratio (S:L) and extraction time (ET). The ranges of their respective investigations are: OT = 40 - 50 ^oC, S:L = 1:20 - 1:60 and ET = 45 - 180 min. Likewise, the responses considered were the extract yield (EY), total phenolic content (TPC) and antioxidant activity (AA). The selected process parameters and investigated responses are consistent with similar works on bioactive antioxidant extractions from plant sources [4,7]. Table 3 shows the different experimental runs (OT, S:L and ET) prescribed by the BBD-RSM and their corresponding values of EY, TPC and AA.

It should be noted that the 3³ BBD produced 17 experimental runs. The range of values of EY = 9.28 – 21.39%, TPC = 18.58 - 69.35 mg GAE/g d.w and AA = 3.51 - 3.77 µM AAE/g. The experimental run 11 (OT = 40 ^oC, S:L = 1:60 and ET = 45 min), run 2 (OT = 30 ^oC, S:L = 1:20 and ET = 112.50 min) and run 4 (OT = 30 ^oC, S:L = 1:60 and ET = 112.50 min) possessed the least EY, TPC and AA respectively. Conversely, the highest EY and TPC and AA belonged to run 5 (OT = 40 ^oC, S:L = 1:60 and ET = 180 min), run 11 (OT = 40 ^oC, S:L = 1:60 and ET = 45 min) respectively. Experimental runs 6, 7 and 16 possessed AA of 3.67 µM AAE/g. The reported range of EY, TPC and AA are consistent with previous studies on HAE-T[7].

3.3. BBD-RSM modelling, model adequacies and statistical analysis

The BBD-RSM was employed to describe the relationship between the studied process parameters and the responses investigated. Hence, predictive mathematical models were developed for EY, TPC and AA as function of OT, S:L and ET with the assistance of Design Expert software. The coded form of the BBD-RSM quadratic predictive equations that relate the TPC, EY and AA with the operating parameters of OT (denoted as A in the equations), S:L (denoted as B) and ET (denoted as C) are presented in Eq. (9), Eq. (10) and Eq. (11) respectively.

$$\begin{array}{l}\mathrm{T}\mathrm{P}\mathrm{C}=+39.16-0.66\mathrm{A}-1.18\mathrm{B}+22.92\mathrm{C}-0.90{\mathrm{A}}^2+2.33{\mathrm{B}}^2+4.15{\mathrm{C}}^2-1.72\mathrm{A}\mathrm{B}-\\ 1.45\mathrm{A}\mathrm{C}+0.38\mathrm{B}\mathrm{C}\end{array}$$

(9)

$$\begin{array}{l}\mathrm{E}\mathrm{Y}=+15.94-2.20\mathrm{A}+4.08\mathrm{B}+1.13\mathrm{C}-1.29{\mathrm{A}}^2+0.047{\mathrm{B}}^2-1.25{\mathrm{C}}^2-1.03\mathrm{A}\mathrm{B}-\\ 0.75\mathrm{A}\mathrm{C}+2.38\mathrm{B}\mathrm{C}\end{array}$$

(10)

$$\begin{array}{l}\mathrm{A}\mathrm{A}=+3.70+0.055\mathrm{A}+7.25\mathrm{x}{10}^{-3}\mathrm{B}-0.055\mathrm{C}-0.040{\mathrm{A}}^2+0.030{\mathrm{B}}^2-\\ 9.85\mathrm{x}{10}^{-3}{\mathrm{C}}^2+9.75\mathrm{x}{10}^{-3}\mathrm{A}\mathrm{B}+0.025\mathrm{A}\mathrm{C}-2.5\mathrm{x}{10}^{-4}\mathrm{B}\mathrm{C}\end{array}$$

(11)

Table 4 contains the ANOVA results for the developed predictive quadratic equations. The models’ coefficients and parameters for assessing their adequacies are as well presented in the table. Table 3 shows that the F-values of the developed EY (32.58), TPC (4947.4) and AA (22.86) models are significant (p values < 0.05). These models also possessed non-significant (p values >0.05) lack of fit indicating they are well-developed and are all capable of predicting the observed experimental data with high accuracy. The R² values for the developed models are also appreciably high (EY = 0.9998. TPC = 0.9767 and AA = 0.9671) and hence indicate high effectiveness and capability of the models in describing the laboratory data. The Pred. R² values of all the models (EY = 0.7077, TPC = 0.9994 and AA = 0.8009) were in close agreement with their respective Adj R² values (EY = 0.9467, TPC =0.9996 and AA =0.9248).

Table 4 further shows the contributions of each term in the equation to the overall predictability of the developed models. The negative (-) and positive (+) sign in the model equations indicated a decrease and increase contributions respectively. Therefore, in the EY BBD-RSM predictive model, the linear terms of OT, S:L and ET are all significant (p < 0.05) and contributed negatively, positively and positively respectively to the overall predictability of the model. The quadratic term of OT (OT X OT) and ET (ET X ET) are significant and they both caused a reduction in the EY model. The quadratic effect of S:L resulted to an insignificant (p > 0.05) increase in EY model. The interactive effect of OT and S:L (OT X S:L) was insignificantly (p > 0.05) negative, however both OT X ET (AC) and S:L X ET (BC) produced a positive contribution in EY model, although the OT X ET effect was not significant (p > 0.05). Regarding the TPC and AA BBD-RSM models, the linear effect of OT, S:L and ET were positive, negative, negative and positive, positive, negative in the TPC and AA models respectively. The linear effects of OT, S:L and ET were significant (p < 0.05) in the TPC model, however, only the linear effects of OT and ET contributed significantly in the AA model. The OT X OT, S:L X S:L and ET X ET significantly (p < 0.05) caused reduction, increment and increment respectively in the developed TPC model. However, the ET X ET was the only insignificant (p > 0.05) quadratic effect in the AA BBD-RSM model. The OT X OT and S:L X S:L effects were significantly (p < 0.05) negative and positive respectively in the AA model. All the interactive effects of OT X S:L, OT X ET and S:L X ET were significant in the TPC model and the respectively produced a negative, negative and positive contributions in the model. Among all the interactive effects, only the OT X ET was significantly (p < 0.05) positive, other interactions such OT X S:L and S:L X ET were insignificant (p > 0.05) negative in the AA BBD-RSM model.

The prediction of the observed experimental data by the developed TPC, EY and AA BBD-RSM models are presented in Figure 2 (a), (b) and (c) respectively. The parity graphs showed excellent prediction of the observed data since the observed experimental and model predicted data fell close to the diagonal line [7].

3.4. Data statistics and multi-gene genetic programming modelling

The properties and variability of both the input and output data were determined prior to MGGP modelling via descriptive statistical analysis. The input and output data were the heat-assisted extraction operating parameters and response parameters respectively as presented in Table 3. Therefore, the input parameters were the OT (^oC), S:L (g/mL) and ET (min) while the output response parameters were the EY (w/w %), TPC (mg GAE/g d.w) and AA (µM AAE/g). The relevant data characteristics that were assessed include data mean, standard error, median, mode, and standard deviation. Others were sample variance, data kurtosis, skewness, range, minimum, maximum, and sum. Seventeen (17) data points of each process parameters (OT, S:L and ET) and response parameters (EY, TPC and AA) (which totaled 102 data population) were used for the construction of the MGGP-based models. Table 5 showed the summary of the descriptive statistics of the data used for the MGGP modelling.

The mean of the data set ranged from 3.6905882 to 41.7864705 with the minimum and maximum belonging to AA and TPC respectively. Also, the mean value of OT, S: L, ET and EY are 40^oC, 40 g/mL, 112.5 min and 14.7682352% respectively. Similarly, the standard error, median, mode, and standard deviation are in the range of 0.0153745 - 11.5761544, 3.69 - 112.5, 3.67 - 112.5, 0.0633907 - 47.7297077 and 0.0040183 - 47.7297077 respectively. The sample variance of both the input and output parameters was between 0.0040183 and 2278.1250. The sample variance of 0.0040183 belonged to AA and it indicated minimal variation while variance of 2278.1250 for ET indicated wide variation. The variance of 271.2990492, 200, and 50 for TPC, S:L and OT also implied wide data variations. The data kurtosis, a statistical parameter that measures the peaked or flatness of a data distribution, was in the range of -0.4367377 - 3.1818297. The kurtosis of - 0.4367377, - 0.6261143, - 0.7428571, - 0.7428571, and -0.7428571 for EY, TPC, OT, S:L and ET respectively, indicated that the data were normally distributed (kurtosis value < +1)) while a value of +3.1818297 for AA implied that the sample population is peaked (kurtosis value > +1) [4]. Similarly, the data skewness (a measure of symmetry of data distribution) was in the range of -1.2675914 - 0.4519267 for all sample population. The skewness coefficients of zero (0) for OT, S:L and ET indicated that the data distributions for these variables were not skewed while a coefficient of -1.2675914 for AA indicated a negatively skewed data distribution. Numerical values of other statistical parameters such as the range, minimum, maximum and sum of the input and output data populations are also indicated in Table 5 and are in the range of 0.26 - 135, 3.51 - 45, 3.77 - 180, and 62.74 - 1912.5 respectively.

The MGGP model structures were optimized for the prediction of TPC, EY and AA as a function of processes extraction variables of operating temperature (OT), solid to liquid ratio (S:L) and extraction time (ET). The MGGP model structure optimization study was purposely carried out to achieve optimum model parameter settings that enable robust learning of the heat-assisted extraction process data (both input and output data presented in Table 3) for finest prediction. Therefore, the population size (PS) and number of generation (NG), which were the two significant parameters that determine the model structural predictive effectiveness to a very large extent [9], were each investigated in the range of 100 – 500. The reported ranges for the parameters investigation were determined based on several attempts of simulation experiments in the preliminary studies. Other model settings are as displayed in Table 1. Table 6 presents the summary of the results of the MGGP modeling optimization studies for the prediction of responses TPC, EY and AA. The R² values for the prediction of TPC, EY and AA were in the range of 0.9858 - 0.9998, 0.9157 - 0.9936 and 0.7595 - 0.9622, respectively.

The parameter settings combination for the optimum prediction of TPC, EY and AA are population size and number of generation of 500 and 250 (R² = 0.9998), 250 and 500 (R² = 0.9936) and 250 and 250 (0.9622), respectively. Therefore, these determined optimum parameter settings for the MGGP model prediction of TPC, EY and AA were henceforth used for the simulation.

Figure 3 presents the graphs (in the form of log Fitness vs. number of generation) of the training process for the best MGGP model structure of TPC (PS = 500; NG = 250), EY (PS = 500; NG = 250) and AA (PS = 250; NG = 250).

It is clear from the figure that the best fitness for the training of the models occurred at approximately 241^th, 498^th, and 245^th generation for the TPC, EY and AA simulations respectively. The fitness function assesses the evolved expressions to determine the most optimal encoded expressions [35] and was obtained by minimizing the root mean square error (RMSE). The prediction errors (also refers to as the fitness) for TPC, EY and AA models were found to reduce as a function of number of generation (NG) and were eventually stable at the optimum NG for each model. The decreases observed in prediction error as training progresses is an indication that the MGGP algorithm learnt from the training data rather than memorizing it and hence no data over-fitting occurred in the course of algorithm training [4]. Also, the ability of all the models to regain permanently decrease in fitness indicated that respective MGGP was able to overcome local minimum and adequately settled in global minimum during training process simulations [4]. The best prediction error for the TPC, EY and AA models were 0.0202, 0.0288, and 0.0011, respectively.

The set of possible solutions capable of predicting the process responses of TPC, EY and AA as a function of set of process variables of OT, S:L and ET are presented in Figure 4 (a), (b) and (c) respectively. Figure 4 showed that Pareto fronts (the circles with green coloration) existed for TPC, EY and AA responses.

The Pareto front represents solutions that outperform all other solutions in both model effectiveness and complexity (measured by the number of nodes in the genetic programming tree) simultaneously. Figure 4 revealed that there are 13, 15 and 14 Pareto solutions for the TPC, EY and AA responses, respectively. However, the green circles with red scribes in Figure 4a, b, and c represent the best Pareto front solutions for the TPC, EY and AA responses respectively, since they possessed the required minimum prediction error which is an indication of high effectiveness in process responses prediction. The mathematical models (best Pareto front solutions) relating the TPC, EY and AA responses to the investigated process variables of OT, S:L and ET are represented by Eq. (12), (13) and (14), respectively. All of the MGGP equations are non-linear in nature and each consisted of five genes and one bias.

$$\begin{array}{l}\mathrm{T}\mathrm{P}\mathrm{C}=7.048\ast {10}^{-4}\left({\mathrm{B}}^2-8.811\mathrm{C}\left(\mathrm{A}-5.1\right)\right)\left(9.149\mathrm{A}+\mathrm{B}-32.2\right)-6.21\ast {10}^{-3}\mathrm{B}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-7.657\ast {10}^{-2}\mathrm{A}+5.215\ast {10}^3\mathrm{C}\left(15.51\mathrm{C}-1.455\right)-9.833\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\ast {10}^{-3}\mathrm{A}\left(9.027\mathrm{A}{\mathrm{C}}^2+33.11\mathrm{A}\mathrm{C}+\mathrm{B}\right)+4.509\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\ast {10}^{-3}\mathrm{C}\left(\mathrm{A}+1.8\mathrm{B}-\mathrm{C}\left(9.149\mathrm{A}+\mathrm{B}\right)-14.52\right)\left(-{\mathrm{C}}^2+\mathrm{A}\mathrm{C}+5.1\mathrm{A}-\mathrm{B}\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+7.048\ast {10}^{-3}\mathrm{A}\mathrm{B}\mathrm{C}-11.17\mathrm{C}\left(8.811\mathrm{C}-3.876\right)\left(\mathrm{A}-12.55\mathrm{A}{\mathrm{C}}^2\right)+163.7\end{array}$$

(12)

$$\begin{array}{l}\mathrm{E}\mathrm{Y}=12.09\mathrm{C}\left(\mathrm{B}\mathrm{C}-\mathrm{B}+\mathrm{C}\left(\mathrm{B}-9.6870\right)\right)-12.09\mathrm{C}-1.633\\ \ast {10}^{-4}{\mathrm{A}}^2\left(\mathrm{B}+\mathrm{C}-9.6870\right)-2.198\ast {10}^{-2}\\ +1.616\mathrm{C}\left(\mathrm{A}+5\mathrm{B}+\mathrm{C}-2.3823\right)-1.099\ast {10}^{-2}{\mathrm{A}}^2\\ +1.616{\mathrm{C}}^2\left(\mathrm{A}+4\mathrm{B}\right)1.099\ast {10}^{-2}\mathrm{A}\left(\mathrm{B}-2.382\right)\\ -\frac{4.127{\mathrm{A}}^2\mathrm{C}\left(\mathrm{C}-4\mathrm{B}+\mathrm{A}\left(\mathrm{A}+5.5470\right)\left(\mathrm{C}+\mathrm{A}\mathrm{C}-3.5975\right)\right)}{{10}^5}-6.532\\ \ast {10}^{-4}\mathrm{B}\mathrm{C}\left(\mathrm{A}+4\mathrm{B}+5\mathrm{B}\mathrm{C}\right)-1.099\ast {10}^{-2}\mathrm{B}\mathrm{C}\left(\mathrm{C}-\mathrm{A}+2\mathrm{B}\mathrm{C}\right)+11.22\end{array}$$

(13)

$$\begin{array}{l}\mathrm{A}\mathrm{A}=2.085\ast {10}^{-1}\mathrm{A}+205.6\mathrm{C}-2.815{\mathrm{A}}^2{\mathrm{C}}^2-6.469\ast {{10}^{-2}\mathrm{A}}^2{\mathrm{C}}^2+4.448\ast \\ {10}^{-2}{\mathrm{A}}^3{\mathrm{C}}^2+7.491\ast {10}^{-2}{\mathrm{A}}^3{\mathrm{C}}^3-1.022\ast {10}^{-2}{\mathrm{A}}^4{\mathrm{C}}^3-6.272\mathrm{A}\mathrm{C}+2.942\ast {10}^{-3}\mathrm{B}\mathrm{C}+\\ 30.93\mathrm{A}{\mathrm{C}}^2-2.942\ast {10}^{-3}{\mathrm{A}}^2\mathrm{C}-3.178\mathrm{B}{\mathrm{C}}^2-3.178\mathrm{B}{\mathrm{C}}^2-1.783\ast {10}^{-3}{\mathrm{A}}^2-1679.0{\mathrm{C}}^2-\\ 3.178\ast {10}^{-2}{\mathrm{C}}^3-3.386\end{array}$$

(14)

The significance (p < 0.05) of each gene in the developed MGGP-based mathematical models for the prediction of TPC, EY and AA values can be visualized in Figure 5 (p-value vs. genes and bias). The structures of the equations, as well as their respective statistical model adequacy parameters are presented alongside in Figure 5.

Figure 5 showed the relative importance (based on the measure of probability values) of each gene and bias that made up the structure of the evolved MGGP equations. All the predictive MGGP-based models for TPC, EY and AA prediction of bioactive antioxidants from P. emblica leaves, possessed interestingly high R² (model R² for TPC, EY and AA are 0.9998, 0.9936 and 0.9622 respectively) and Adj R² (Adj R² for TPC, EY and AA are 0.9997, 0.9907 and 0.9451) values which indicate that the models are a good fit for the observed data and that they explain a significant portion of the variability in the dependent variables. Considering the evolved MGGP equation for the prediction of TPC (Eq. (11)), the probability all the genes and the bias were less than 0.0002 which implied that they are all significant (p < 0.05). However, the most important genes are genes 2, 3, 4 and 5. The bias was also a significant model term for MGGP-based predictive model for TPC as indicated in the figure. Similarly in Figure 5 (b), all the genes and the bias were significant part of MGGP-based predictive model structure for EY prediction since all the probability values of the evolved genes and bias were less than 0.00015. However, of significantly high importance in the evolved EY MGGP-based predictive model are the bias, gene 1, gene 2 and gene 5. Analysis of Figure 5 (c) also showed that the MGGP-based AA model have significant structure (genes and bias) for the prediction of AA values. As indicated on the graph (Figure 5 (c)), the probability values of all genes and bias in the structure are less than 0.000015, however, gene 1 is the most significant gene in the structure.

The ability of the MGGP-based models to predict the laboratory observed TPC, EY and AA data as function of process variables of OT, S:L and ET are presented in Figure 6 (a), (b) and (c), respectively.

It is obvious from Figure 6 that the MGGP-based models were able to predict the observed laboratory data for TPC, EY and AA satisfactorily. The parity graphs showed that both the observed and predicted data were clustered on the diagonal line which is an indication of good data predictive strengths of the MGGP-based models [4]. The models’ RMSE statistical parameter values were minimal (RMSE for TPC, EY and AA MGGP-developed models are 0.0202, 0.0288 and 0.0011) while all the models’ R² values were close to unity which implied perfect prediction of the models [4].

3.5. Comparison of the BBD-RSM and MGGP - based predictive models

The efficiency of RSM and MGGP is compared using their respective statistical indicators and the result is summarized in Table 9.

Adeyi et al. [4] reported that models with higher R² value implied better performances. Here, Table 9 depicts that the two models (i.e. RSM and MGGP model) performed well, however, the R² – value of RSM model is marginally greater than that of MGGP model and means a better model performance. RSM models are statistical models derived from regression analysis to fit a model to the data. It means that RSM models are parametric and therefore requires a priori knowledge of the functional form of the model [36]. However, parametric models are computationally efficient and can be utilized to model complex systems with a small number of data. On the other hand, MGGP is non-parametric and thus not require a priori knowledge of the functional form of the model, however, it is computationally intensive and evolves more accurate models as data utilized increased [37]. Therefore, the slightly improved performance of RSM model to MGGP model can be well related to small data set utilized, hence, RSM modeling approach is recommended for further model implementation.

3.6. Numerical optimization and validation

The desirability algorithm that is present in the Design Expert software was used for achieving the optimization of the process parameter for the recovery of bioactive antioxidants from P. emblica leaves. The optimization procedure is in accordance with the work of Adeyi et al. [16]. The goal of the optimization scheme was to determine a set of process parameters that maximized the TPC, EY and AA of the bioactive extract. Table 5 presents the selected goal, weight, and importance for both the process parameters and responses. The table further shows the numerical range of search for the global optimum parameters of processing bioactive antioxidants recovery from P. emblica leaves.

Table 6 summarized all the thirteen (13) solutions presented based on the combined desirability values. The scale of the desirability is between 0 and 1, with desirability value of 1 adjudged the best [4].

The desirability values of the all the presented solutions ranged from 0.766 to 0.854. These solutions were ranked in ascending order of preference by the software with the first solution on the list, the best and selected. The 13^th solution on the table was the least preferred, although the solution also possessed appreciably high desirability value (0.766), the value was the least among the solutions provided by the algorithm. Hence, the first solution on the list with desirability value of 0.854 was selected as the best optimum solution. This desirability value of 0.854 compares well with other desirability values of selected optimum solutions in the literature. For instance this present desirability value is higher that than the desirability value of selected optimum solution in the work of Adeyi et al., [4] during the HAE process optimization of bioactive extract recovery from Carica papaya leaves.

Figure 3 is the optimization ramp for the best suggested solution (first) with combined desirability value of 0.854. The figure shows that the optimization search was within the range of investigated process parameters and observed laboratory response data. Therefore, the process parameters that simultaneously gave the optimum EY of 21.6565%, TPC of 67.116 mg GAE/g and AA of 3.68583 µM AAE/g were OT of 41.61 ^oC, S:L of 1:60 g/mL and ET of 180 min.

The validation experiment was conducted in the laboratory to ascertain the selected predicted optimum in the laboratory. Therefore, 1 g of P. emblica leaves was mixed with 60 mL of 50% ethanol-water mixture in a beaker and heated to approximately 42 ^oC by utilizing a water bath for the duration of 180 min. After the completion of the experiment, the leave fibers were separated from the extract through centrifugation and the EY, TPC and AA were quantified according to Section 2.7, 2.8 and 2.9 respectively. The results obtained for the validation experiment were EY = 22.31 %, TPC = 69.612 mg GAE/g and 3.72 µM AAE/g. The relative standard deviation (RSD) was computed to compare the experimental validated result with the predicted response result. It was found that the RSD between the validated and predicted values of EY, TPC and AA were 2.67%, 7.45% and 4.68% respectively. This indicated that these values are similar since the RSD between the validated and predicted response values are less than 10 [16]. This similarity in the validated and predicted results showed that the developed BBD-RSM models for EY, TPC and AA are effective, well-fitted, robust and capable of predicting the process of bioactive extract recovery from P. emblica leaves.

3.7. Phenolic fingerprints of P. emblica leaf extract

HPLC profiling of phenolic compounds in plant extracts is of paramount importance for the identification, and characterization of these bioactive compounds. It provides valuable insights into the chemical composition, quality, and therapeutic potential of plant extracts. Additionally, HPLC profiling can aid in the evaluation of the stability and degradation kinetics of phenolic compounds under different conditions, ensuring the preservation of their bioactivity. Therefore, the HPLC profiling of phenolic compounds in P. emblica leaf extract was for the purpose of identifying the phenolic compounds with potential therapeutic effects for the basis of future industrial production and techno-economic analysis of the P. emblica leaf extract. Hence to this end, eight (8) phenolic compounds with established bioactivities were used as standards for the HPLC profiling of the P. emblica leaf extract. These compounds were compared with the content of the extract using the similarities in retention factor (RF). Figure 4 is the HPLC fingerprints of P. emblica leaf extract. As shown, Figure 4 is characterized by different phenolic compounds and wide disparities in their corresponding areas. The figure revealed six (6) phenolic compounds that have comparable RT with the HPLC phenolic standards used as the baseline of identification. The identified phenolic compounds in the P. emblica leaf extract with corresponding RT are betulinic acid (RT = 2.439 min), gallic acid (RT = 3.063 min), chlorogenic acid (RT = 3.541 min), caffeic acid (RT = 4.055 min), ellagic acid (RT = 5.825 min), and ferulic acid (RT = 7.684 min).

The bioactivities of the identified phenolic compounds in the extract were interesting and pointing to the overall therapeutic potential of P. emblica leaf extract. For instance, betulinic acid exhibits a wide range of therapeutic properties and studies have shown that it induces apoptosis (programmed cell death) in cancer cells, making it a promising candidate for cancer treatment [38]. Gallic acid possesses numerous health benefits and recent investigations showed that it has potentials to inhibit tumor growth, reduce oxidative stress, and exert protective effects against various diseases [39]. Chlorogenic acid has been studied for its potential in preventing skin tumorigenesis, modulating MAPK and NF-κB pathways, and ameliorating oxidative stress [40] while caffeic acid has been investigated for its potential in inhibiting atopic dermatitis-like skin inflammation and synergistic antioxidant activity when combined with other phenolic acids [39]. Both ellagic acid and ferulic acid have been shown to exhibit antioxidant, anti-inflammatory, anticancer, antimicrobial, and antimutagenic effects [41,42].

3.8. Reliability Assessment of BBD-RSM based Predictive Models

Model reliability assessment is the process of evaluating and determining the trustworthiness and accuracy of a model’s predictions or outputs. One of the various available techniques and methodologies to gauge the model’s performance and identify its strengths and limitations is via Monte Carlo simulation. The main purpose of Monte Carlo simulation in the context of model reliability assessment is to quantify uncertainty, validate model performance, and estimate the potential range of outcomes. Another key aspect of Monte Carlo simulation is its ability to conduct robust sensitivity assessment where various inputs are systematically varied to determine their impact on the model’s outputs which helps in identifying which input parameters have the most significant influence on the results and helps in understanding the robustness of the model. In the present investigation, the developed BBD-RSM based predictive models were used for this analysis due to their relative superiority in predicting the heat-assisted extraction process outputs of TPC, EY and AA, relative to the MGGP-based models as explained in Section 3.5. Figure 7 (a), (b) and (c) show the split views of the probability distributions, cumulative frequency and reverse cumulative frequency curves of the analyses conducted in Oracle Crystal ball software for the prediction of experimental observed data (in Table 3) of TPC, EY and AA respectively, as function of input variables of OT, S:L and ET in order to ascertain the robustness of the constructed BBD-RSM based predictive models. The outcomes’ data statistics and the best fit model for the distribution were also incorporated in the presented split view figures and positioned at the upper right and down respectively. Figure 7 shows that the BBD-RSM models for TPC, EY and AA are capable of predicting the respective outcomes as function of OT, S:L and ET within the range of their respective experimental outcomes. The data outcomes kurtosis and skewness (statistical measures that provide information about the shape and distribution of a dataset) for TPC, EY and AA are 2.77 and 0.2508; 3.64 and 0.5979; and 5.40 and -1.14, respectively. The distribution curves for the TPC, EY and AA were excellently modeled by Weibull (Anderson-Darling coefficient = 308.1560), Lognormal (Anderson-Darling coefficient = 38.1254) and Min Extreme (Anderson-Darling coefficient = 34.7595), respectively. As indicated in Figure 7, the percentage certainty of the developed BBD-RSM models in predicting the observed experimental data range of 18.58 – 69.35 mg GAE/g, 9.28 – 22.14% and 3.51 - 3.77 µM AAE/g (in Table 3) for TPC, EY and AA, is 99.985%, 97.569% and 98.661%, respectively. The high BBD-RSM models’ individual prediction certainty is an indication of their respective high reliability and robustness.

These values compare well with the prediction certainty of the Sierrathrissaleonensis cracker drying effective moisture diffusivity D-Optimal-RSM predictive model (99.831%) [14], biodiesel yield CCD-RSM predictive model (73.509%) (Oke et al., 2022) and techno-economic MGGP predictive models (MGGP-CAnysP APR = 99.980% and MGGP-CAnysP UPC = 98.477%) [14]. The dynamic sensitivity charts, which identify input parameters that have the most significant influence on the predictions of the developed BBD-RSM models, are presented in Figure 8. Here, the contribution of process variables of OT, S:L and ET to the variance in the prediction of BBD-RSM models for TPC, EY and AA were assessed. Figure 8 shows that process variables contributed differently to the perturbation in the developed BBD-RSM models predictions. In the dynamic sensitivity graphs, all the bars to the right hand side indicate positive contributions (increase in response variable value with an increase in process variable value) while the bars to the left hand side signify negative contributions (increase in process variable value with decrease in response variable value).

Also the length (measured by percentage) of a sensitivity bar determines its magnitude and relative importance. Hence, a long bar (to either side) is relatively significant than correspondingly shorter bar, with 0% assigned to a no-significant effect. A high significant parameter should be better controlled (or more accurately measured) in order to improve the model’s predictability. Therefore, analysis of Figure 8 (a) showed that S:L was the process variable with the most significant importance to the variance in TPC BBD-RSM model prediction. The ET and OT did not seem to influence profoundly the predictability of TPC BBD-RSM model. In the same vain, Figure 8 (b) showed that the ET and S:L were the process variables with highest and least significance, respectively to the variance in the EY BBD-RSM model prediction. Also the order of process variable significance (Figure 8 (c)) to the perturbation in the predictability of AA BBD-RSM model is S:L > OT > ET. In numerical terms however, S:L had positive contribution (+99.7%), ET had negative contribution (-0.28%) while OT did not contribute significantly (+ 0.02%) to the perturbation in the predictability of TPC BBD-RSM. Likewise, ET contributed positively (+76.2%), OT contributed negatively (-19.8%) while S:L contributed positively (+4%) to the variance in EY data prediction by the EY BBD-RSM model. The S:L, OT and ET contributed -52%, +46.7 and +1.3% respectively to the variance in AA value predictability by the constructed AA BBD-RSM model.

4.0. Conclusion

It is feasible to recover bioactive antioxidants from Phyllantus emblica leaves using the heat-assisted extraction technology and water-ethanol solvent mixture. The best ethanol-water mixture ratio for optimal recovery of bioactive extract was 50:50. All the investigated heat-assisted extraction process parameters had notable impacts on the bioactive extract’s total phenolic compound (TPC), yield (EY) and antioxidant activity (AA). The process parameters that simultaneously gave the optimum EY of 21.6565%, TPC of 67.116 mg GAE/g and AA of 3.68583 µM AAE/g were OT of 41.61 ^oC, S:L of 1:60 g/mL and ET of 180 min. The HPLC extract characterization revealed that the bioactive constituents of the recovered extract were betulinic acid, gallic acid, chlorogenic acid, caffeic acid, ellagic acid, and ferulic acid. Both response surface methodology (RSM) and multi gene genetic programing algorithm (MGGP) were able to model the relationships between the investigated heat-assisted extraction process parameters and observed process responses of TPC, EY and AA with high satisfaction. The developed RSM models (BBD-RSM) were however, marginally better than their corresponding MGGP based models for the prediction of TPC, EY and AA data. The BBD-RSM based models were highly reliable and were able to predict the laboratory observed data of TPC, EY and AA as function process variables with high certainty. The certainty values for the developed BBD-RSM models to predict observed TPC, EY and AA values were 99.985%, 97.569% and 98.661%, respectively. The process parameter of S:L had strong positive contribution (+99.7%) while OT, and ET did not seem to have strong effect (+ 0.02% for OT and -0.28% for ET) on the TPC values prediction of the BBD-RSM model. In the same manner, ET contributed positively (+76.2%), OT contributed negatively (-19.8%) while S:L contributed positively (+4%) to the variance in EY data prediction by the EY BBD-RSM model. The order of contribution of the process variables to the variance in AA value prediction by the constructed AA BBD-RSM model was S:L > OT > ET.