Martin M. Kasina¹, Koske Joseph², Mutiso John^2
1Department of Mathematics and Statistics, Machakos University, Machakos, Kenya.
²Department of Mathematics, Physics and Computing, Moi University, Eldoret, Kenya.
DOI: 10.4236/ojop.2020.93005   PDF    HTML   XML   498 Downloads   2,253 Views   Citations

Abstract

Despite the increased recognition of the nutritional value of the Oyster mushroom, coupled with its ability to tolerate a wide range of climatic conditions, its production is still at infancy stage with low adoption rate in Kenya. The low uptake could be attributed to the cost of spawns or lack of skills for spawns preparations coupled with poor knowledge on oyster mushroom consumption benefits. The objective of this study was to optimize Pleurotus ostreatus (Oyster mushroom) spawns production. To achieve the objective, the spawns propagation was optimized by varying the temperature level, sterilization time and culture media concentration in order to establish the feasible levels which minimized the days of mycelium full development using central composite designs. Based on the study findings, 26.29˚C, 17.36 minutes and 60.95 g/L of temperature level, sterilization time and culture media concentration levels respectively minimized the days to full coverage of mycelium in a petri dish. Central composite designs for controlling temperature, sterilization time and culture media concentration were recommended for spawns maximum production. A further research on multiple response optimizations aimed at achieving resistance to bacterial diseases and yield by varying the strain in the culture were recommended.

Keywords

Spawns, Mycelium, Colonize, Optimal Levels, Second Order Model, Response surface

Share and Cite:

1. Introduction

Central Composite Design (CCD) is a common form of response surface methodology. The designs are mostly used when the design appropriately fits sequential experimentation since the designs obtain information from a correctly planned factorial experiment which can accommodate up to five levels per factor [1].

The CCD is composed of three design points: square or edge points in two level designs (±1), star points at $\pm \alpha$ ; $\left|\alpha \right|\ge \text{1}$ that take care of the quadratic effect and then the third one is the centre points [2].

The square or cube points consist of 2^k factorial design (±1, ±1, …, ±1), where k is the number of independent variables in which the design can be replicated $n_f$ times. The star points consist of 2k units on the axis of each factor at a distance α, from the centre of the design [(±α, 0, …, 0), (0, ±α, 0, …, 0)], whose selection is based on the orthogonality and rotatability criterion, and takes one observation at each of the vector $\pm \alpha e_i$ can be replicated $n_s$ times, where $e_i$ is the i-th Euclidean unit vector and $\alpha >0$. The centre points designated as “0” points and can be replicated $n_c$ times. Thus by letting n denote the total number of experimental runs in the CCD, based on k design factors, the runs size can be determined by using Equation (1)

$n=2^{k-p}{n}_f+2k{n}_s+n_c$ (1)

where:

n is the desired sample size.

$2^{k-p}{n}_f$ is factorial component of the design; $n_f$ is the number of replication times; k is the number of factors varied and p is the number of factors subtracted from k.

$2k{n}_s$ is the axial component of the design replicated $n_s$ times.

$n_c$ is the number of times centre points that were replicated.

There are various types of central composite designs depending on where the star points are placed, namely: circumscribed, inscribed and face centered central composite designs. Circumscribed Composite Designs (CCC) has a circular, spherical or hyperspherical symmetry. In the circumstances where the limits specified for settings are truly the limits, then it is referred to as the Inscribed Composite Design (CCI) which uses the factor settings as the star points and creates a factorial or fractional factorial design within those limits. Then lastly Face Centered Composite Design (CCF) has an alpha of 1, such that the axial points are at the center of each face of the factorial space, hence, $\alpha \pm 1$. Figure 1 shows the graphical representation of the three variants of central composite designs;

circumscribed (CCC), face centred (CCF) and inscribed (CCI) with two level factorial experiments for three factors respectively.

1.1. Spawns and Oyster Mushroom Production

Mushrooms are macrofungi with a distinctive fruiting body which can be either epigeous or hypogeous. The macrofungus has fruiting body large enough to be seen with the naked eye and can be picked up by hand [3].

Spawn, which is the reproducing part of the mushrooms, should be multiplied in a pure culture media in a sterilized environment. The pure culture can either be raised by tissue culture or spore culture. In tissue culture a well grown mushroom with membrane covering the gills is selected, from which a small bit of mushroom from gill portion is taken using forceps and inoculated on Potato Dextrose Agar (PDA) or Malt Extract Agar (MEA). Under spore culture method, the spores are collected from well-developed fruiting body by “spore mapping technique” and then the spores are inoculated to the PDA or MEA slants as in tissue culture under aseptic condition. The mycelia (spawns) are supposed to cover the entire surface after which they are placed in a sterilized jars or bottles with a boiled grain seeds like wheat or millet for further multiplication. It is at this stage when the seeds can be used as mushroom spawn and can produce mushrooms if mixed with a well compost substrate.

1.2. Problem Statement

Despite the increased recognition of the nutritional value of the Oyster mushroom, coupled with its ability to tolerate a wide range of climatic conditions, its production in Kenya is still at infancy stage with low adoption rate [4]. The low uptake could be attributed to the cost of buying spawns which could be substituted by the acquisition of the skills to prepare them, coupled with poor knowledge on its production and consumption benefits.

Past studies indicates that, beside the correctly compost substrate, the oyster mushroom does well under certain levels of temperature, humidity, light and ventilation, while the quality of spawn is affected by the mother culture [5], media culture concentrations and spawn substrate preparations [6]. Other several studies indicate that oyster mushrooms yield depends on the type of substrate used and the quality of the spawns [7] [8].

Hence this study focused on determining the optimal level of temperature, sterilization time for media culture and the media culture concentration settings, that minimized the time (days) to form mycelium culture establishment during the oyster mushroom spawn preparation, by using central composite design. This would reduce the process, production costs and narrow the gap of opportunistic competitors’ invasion, hence improve on the oyster mushroom’s spawn quality, production and consumption levels.

2. Literature Review

Central Composite Design has widely been used in various fields of research; it was used in optimizing machine parameters in drilling hybrid metal matrix composites, and it was found that the predicted values and measured values were fairly close, indicating that the developed models could be effectively be used to predict the responses in the drilling of hybrid metal matrix composites, [9]. CCD’s are very popular 2^nd-order designs because; they are extremely simple to use, and allows estimation of all the parameters in a full second-order model coupled with the fact that among the exact (integer) designs, they often have high efficiencies under the commonly used A-, D- and E-optimality criteria [10].

Composite Design and Response Surface Methodology were applied to optimize the conditions of submerged cultivation of the Fomesfomentarius mushroom. Concentration of dextrose, yeast extract and time of cultivation in days were selected and their correlative effect on mushroom multiplication established. The maximum yield of dry weight biomass was 23.74 g/L with 0.8 - 7.5 g/L concentration after 11 days of submerged cultivation [11].

CCD was used to show that, the evaporative cooling of the air surrounding the pileus created convective airflows capable of carrying spores at speeds of centimetres per second. That work revealed how mushrooms tolerate and even benefit from crowding and explained their high water needs. It was evident that spores continuously flow out from thin gaps, even in the absence of external winds [12].

CCD was used to optimize the Degradation of Tylosin from Aqueous solution by Photo-Fenton Reaction. The interaction effects and optimal parameters were obtained and the significance of the independent variables and their interactions was tested by means of analysis of variance (ANOVA) at 5% significance level. The results showed that the concentration of the ferrous ion and pH were the main parameters affecting Total Organic Carbon (TOC) removal, while peroxide concentration had a slight effect on the reaction. The optimum operating conditions to achieve maximum TOC removal were determined. There was a good agreement between the model prediction and experimental results confirming the soundness of the developed model [13].

Response surface methodology was applied to optimize the main factors which significantly affected exopolysaccharide (EPS) production. The concentrations of glucose and peptone were found to be the main effective factors for EPS production by the fractional factorial design (FFD) and central composite design experimental analysis. Verification experiment confirmed the validity with the actual EPS yield as 1.97 g/L, which was 6.29-fold in comparison with the (2.22 g/L) in the original basal medium [14].

Central composite design and response surface methodology were used in the advanced treatment of olive oil processing wastewater using Fenton’s peroxidation, in a second-order polynomial multiple regression model. Analysis of variance (ANOVA) showed a high coefficient of determination (R²) value of 0.902 - 0.998, thus ensuring a satisfactory, adjustment of the second-order regression model with the experimental data [15].

Under this study, the central composite design was used for theoretical examination of optimal temperature, sterilization time and the media culture concentration experimental space that minimized the time (days) to full spawn growth.

3. Methodology

The growth process of the spawns was a function of temperature level, sterilization time and the culture media concentration level as shown in Equation (2)

$y=f\left(x_1,x_2,x_3\right)+\epsilon$ (2)

$x_1$ = Temperature level

$x_2$ = Sterilization Time

$x_3$ = Culture media Concentration level

The objective function was to determine the level of the regressors’ experimental setting that minimized the mycelia full coverage in the petri dish area through the response surface methodologies.

A second order response surface design should be rotatable, to make it possible to extract the most information regarding the dependent variable, and leave the least amount of uncertainty for the prediction of future values [16].

Rotatable design provides the preferred property of constant prediction variance at all points that are equidistant from the design center, hence improving the quality of the prediction. It is a desirable property, especially when there is a need to optimize $\hat{y}\left(x\right)$ over the region of interest. A design is said to be orthogonal if it can provide independent information about the effects of the various terms in the model.

In order to ensure orthogonality and rotatability of the design, circumscribed composite design was preferred. This ensured that any non-allowable operating conditions at two or more of the extremes of the design region were encompassed.

3.1. Rotatability Conditions

The design matrix for a CCD experiment involving k independent variables of L levels, is a matrix derived from the values corresponding to the three types of experimental runs stated in Equation (2) to form a design with N of $L^N$ treatment combinations matrix, such that

$X=\left[\begin{array}{cccc}\pm 1_{1,1},& \pm 1_{1,2}& \cdots & \pm 1_{1,m}\\ ⋮& ⋮& \ddots & ⋮\\ \alpha & 0& \cdots & 0\\ -\alpha & 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & \alpha \\ 0& 0& \cdots & -\alpha \\ 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0_{N,1}& 0_{N,2}& \cdots & 0_{N,m}\end{array}\right]$

usually referred to as the design matrix. Where, the treatment combinations are the points of the design, such that $x_i$ denotes the i^th factor at j^th treatment in a given level. Such a design will be a rotatable design of order d if a response polynomial surface is of Equation (3) form [16].

$\begin{array}{c}{y}_u={\beta }_{0u}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\beta }_{iu}{x}_{ui}}+{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\beta }_{iiu}}{x}_{ui}^2+{\displaystyle \underset{i=1}{\overset{n-1}{\sum }}{\displaystyle \underset{i<j=2}{\overset{n}{\sum }}{\beta }_{iju}}}{x}_{ui}{x}_{uj}\\ +{\displaystyle \underset{i=1}{\overset{n-2}{\sum }}{\displaystyle \underset{i<j=2}{\overset{n-1}{\sum }}{\displaystyle \underset{i<j<k=3}{\overset{n}{\sum }}{\beta }_{ijku}{x}_{iu}{x}_{ju}}}}{x}_{ku}+\cdots +{\epsilon }_u\end{array}$ (3)

which can be so fitted such that the variance, $Var\left[\hat{y}\left(x\right)\right]$ is a function of $r_u^2=x_{1u}^2+x_{2u}^2+\cdots +x_{ku}^2$ only, implying the variance of ${\hat{y}}_u$ is a function of only the distance of x from the origin but not of the direction and therefore constant at all points that are equidistant from the design center.

When the model is of the second degree, that is $d=2$, then such constancy of variance is attained if the design points are chosen to satisfy the following conditions:

${\displaystyle \underset{u=1}{\overset{n}{\sum }}{x}_{iu}}=0,{\displaystyle \underset{u=1}{\overset{n}{\sum }}{x}_{ui}}{x}_{uj}=0,{\displaystyle \underset{u=1}{\overset{n}{\sum }}{x}_{iu}{x}_{ju}^2}=0,\text{for}i\ne j=0,1,\cdots ,k$ (i)

${\displaystyle \underset{u=1}{\overset{n}{\sum }}{x}_{ui}^2}=\text{constant}=N{\lambda }_2,\forall i$ (ii)

${\displaystyle \underset{u=1}{\overset{n}{\sum }}{x}_{ui}^4}=\text{constant}=3N{\lambda }_4,\forall i$ (iii)

${\displaystyle \underset{u=1}{\overset{n}{\sum }}{x}_{ui}^2{x}_{uj}^2}=\text{constant}=N{\lambda }_4\forall i\ne j$ (iv)

where ${\lambda }_3$ and ${\lambda }_4$ are constants of any design.

Hence, ${\displaystyle \underset{u=1}{\overset{n}{\sum }}{x}_{iu}^4}=3{\displaystyle \underset{u=1}{\overset{n}{\sum }}{x}_{iu}^2}{x}_{ju}^2$ for $i\ne j=1,2,\cdots ,k$ (4)

where the summation of the relations in Equation (4) is over the design points $u=1,2,\cdots ,N$.

The choice of axial distance $\alpha$ is based on the region of interest and operability. This study considered rotatable central composite design (RCCD) under which the value of $\alpha$ depends on the number of experimental runs in the factorial

portion of the central composite design, such that $\alpha ={\left(f\right)}^{\frac{1}{4}}$, where f is the

square points in a central composite design, making the design rotatable.

Therefore by using Equation (4) and the rotatability conditions, Equation (5) was obtained [16] [17] [18] [19].

${\alpha }^4=\frac{2^{k-p}{n}_c}{n_s}$ (5)

3.2. The First- Order Model

The first order model aids in getting the direction of the relationship between the response and independent variables. It leads the experimenter sequentially and efficiently along the path of steepest ascent/descent. The first order model is a lower-order polynomial. A first order model with k variables takes the form expressed in Equation (6).

$y={\beta }_0+{\beta }_1{x}_1+\cdots +{\beta }_k{x}_k+\epsilon$ (6)

where $x_i$ and ${\beta }_k$ ’s are the design variables and regression coefficients respectively.

When the curvature appears, the second-order model is used to approximate the response values precisely.

3.3. The Second-Order Model

The central composite design was used to explore the region for fitting in the first and second order models in this study. The prevailing region of exploration was 20˚ - 30˚ degree Celsius for temperature, 10 - 20 minutes of sterilization time and 35 - 65 grams of agar per 1000 cubic centimeter of PDA. To simplify the calculations the independent variables were coded (−1, 1) interval such that given ${\xi }_1$, ${\xi }_2$ and ${\xi }_3$ as the natural variables for temperature, time and culture media concentration levels respectively, then the coded variables should be derived at by using the Equation (7), relations [16] [20].

$x_1=\frac{{\xi }_1-25}{5}$, $x_2=\frac{{\xi }_2-15}{5}$ and $x_3=\frac{{\xi }_3-50}{15}$ (7)

To enable in the identification of the best CCD by creating a grid of all combinations of the design choices, the information function for the second-order (quadratic) model had to be rotatable, the $\alpha$ values necessary for orthogonality and rotatability were computed. Ideally, in order to make an unbiased estimate of pure error, the CCD should comprise of three to five centre points [21]. The central composite design is commonly used to fit the second order response surface model of the form expressed in Equation (8)

$y_u={\beta }_0+{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\beta }_i{x}_{ui}}+{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\beta }_{ii}}{x}_{ui}^2+{\displaystyle \underset{i=1}{\overset{k-1}{\sum }}{\displaystyle \underset{i<j=2}{\overset{k}{\sum }}{\beta }_{ij}}}{x}_{ui}{x}_{uj}+{\epsilon }_u$ (8)

where;

$y_u$ is the response obtained from the u^th combination of factors ( $u=1,2,\cdots ,n$ ),

$x_{ui}$ denotes the level of the i^th factor ( $i=1,2,\cdots ,k$ ) in the u^th run ( $u=1,2,\cdots ,n$ ) of the experiment,

${\beta }_i$ ’s are the model coefficients to be determined; ${\beta }_0$ being a constant, ${\beta }_i$ being the i^th linear regression coefficient, ${\beta }_{ii}$ being the i^th quadratic regression coefficient and ${\beta }_{ij}$ being the ${\left(i,j\right)}^{\text{th}}$ interaction coefficient and

${\epsilon }_u$ ’s are the uncorrelated random errors in the u^th observation with a mean zero and a constant variance ${\sigma }^2$.

3.3.1. The Significant Regression Coefficients

To test the significance of the individual regression coefficient ${\beta }_j$, the P-Value and the t test were used. The t test statistics is based on the t distribution and has the form of Equation (9)

$t=\frac{{\hat{\beta }}_j}{s\left({\hat{\beta }}_j\right)}$ (9)

where ${\hat{\beta }}_j$ is the least square estimator of the parameter ${\beta }_j$ ( $j=1,2,\cdots ,q$ ) and $s\left({\hat{\beta }}_j\right)$ is the estimated standard error of ${\hat{\beta }}_j$. The standard error of each parameter ${\hat{\beta }}_j$ is given by the square root of the diagonal elements of the matrix, $var\left({\hat{\beta }}_j\right)$ where

$Var\left({\hat{\beta }}_j\right)={\sigma }^2{\left(X^{\text{T}}X\right)}^{-1}$

and X is the matrix of values of explanatory variables referred to as the design matrix, such that.

$X=\left[\begin{array}{ccccc}1& x_{1,1}& x_{1,2}& \cdots & x_{1,p-1}\\ 1& x_{2,1}& x_{2,2}& \cdots & x_{2,p-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n,1}& x_{n,2}& \cdots & x_{n,p-1}\end{array}\right]$

3.3.2. Location and Characterization of the Stationary Points

The temperature level ( $x_1$ ), Sterilization time ( $x_2$ ) and Culture media concentration level ( $x_3$ ), that minimized the days of Pleurotus ostreatus spawns development were established by carrying out the partial derivatives from the function $f\left(x_1,x_2,x_3\right)=\hat{y}$ as shown in Equation (10)

$\frac{\partial \hat{y}}{\partial x_1}=\frac{\partial \hat{y}}{\partial x_2}=\frac{\partial \hat{y}}{\partial x_3}=0$ (10)

and by letting $b_0,b_i,b_{ii},b_{ij}$ denote the least square estimators of ${\beta }_0,{\beta }_i,{\beta }_{ii},{\beta }_{ij}$ respectively we have the response and the fitted values taking the form of Equation (11)

$\hat{y}=b_0+b_i{x}_i+b_{ii}{x}_i^2+\underset{i<j}{b_{ij}}{x}_{ij}$ (11)

whereby the second order model in matrix notation is as summarized in Equation (12)

$\hat{y}={\hat{\beta }}_0+x^{1}b+x^{1}Bx$ (12)

such that the fitted values in matrix form are;

$X=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_k\end{array}\right]$, $b=\left[\begin{array}{c}{\hat{\beta }}_1\\ {\hat{\beta }}_2\\ ⋮\\ {\hat{\beta }}_k\end{array}\right]$ and $B=\left[\begin{array}{cccc}{\hat{\beta }}_{1,1}& \frac{{\hat{\beta }}_{1,2}}{2}& \cdots & \frac{{\hat{\beta }}_{1,k}}{2}\\ {\hat{\beta }}_{2,1}& {\hat{\beta }}_{2,2}& \cdots & \frac{{\hat{\beta }}_{2,k}}{2}\\ ⋮& ⋮& \ddots & ⋮\\ sym& \frac{{\hat{\beta }}_{k,2}}{2}& \cdots & {\hat{\beta }}_{k,k}\end{array}\right]$

whereby b is a (k × 1) vector of first order regression coefficients and B is a (k × k) symmetric matrix whose main diagonal elements are pure quadratic coefficients ( ${\hat{\beta }}_{ii}$ ) and whose off-diagonal elements are one-half the mixed quadratic coefficients ( ${\hat{\beta }}_{ij},i\ne j$ )

Therefore at the stationary points $\frac{\partial \hat{y}}{\partial x}=b+2Bx=0$ and the stationary points would be $x_s=\frac{1}{2}{B}^{-1}b$ and the predicted response at the stationary points would be determined by Equation (13)

${\hat{Y}}_s={\hat{\beta }}_0+\frac{1}{2}{x}_s^{1}b$ (13)

The second order model normally takes the maximum, minimum, saddle point or stationary ridge. The coordinates and the nature of each stationary point were determined by carrying out the partial derivatives and represented by using the contour maps.

3.3.3. Canonical Analysis

The canonical analysis objective was to determine the nature of the stationary point and the entire response surface. If the signs of the eigenvalues of $\hat{B}$ are all positive then the stationary point is minimum otherwise it is a maximum response surface if all the signs of the eigenvalues of $\hat{B}$ are all negative. A saddle point would exist if there is a mixed sign, which is both positive and negative.

The canonical equation takes the form of Equation (14)

$\stackrel{⌢}{y}={\stackrel{⌢}{y}}_s+{\lambda }_1{w}_1^2+{\lambda }_2{w}_2^2+\cdots +{\lambda }_k{w}_k^2$ (14)

where $\left\{{\lambda }_i\right\}$ are the eigenvalues or the characteristic roots of the matrix B, while $\left\{w_i\right\}$ are the transformed independent variables. The nature of the response surface is determined from the stationary points and the signs of the magnitudes ${\lambda }_i$. If all ${\lambda }_i$ are positive then $x_s$ is a minimum response surface, otherwise $x_s$ is a maximum response surface if all ${\lambda }_i$ are negative, while mixed sign indicates a saddle response surface and if at least one eigenvalue is zero then it is a ridge surface. The response surface is steepest in the $\left\{w_i\right\}$ direction corresponding to the largest absolute eigenvalue. Canonical analysis is very instrumental in systems of maxima and minima determination in many dimensions and, in particular to identify complicated ridge systems, where direct geometric representation is not possible.

4. Results and Discussion

4.1. Optimization of Spawns Propagation

The operating optimal levels of temperature, sterilization time and culture media concentration that minimized time in days of the mycelia full coverage in a petri dish area were determined through central composite designs. The colonised media in a petri dish is displayed in Figure 2.

4.1.1. Central Composite Design

The central composite design was used to explore the region for fitting in the first and second order models. The second order design of 3 factors; temperature level, sterilization time and culture media concentration at 2 levels were chosen

as the process variables for optimization, by investigating their effect on time to full colonization of a substrate in a petro dish for spawn production of oyster mushroom (Pleurotus ostreatus). The CCD comprised a total of 19 experiments: a full factorial design 8 experiments, 6 axial points, and central points replicated five times. In order to measure the effect of incubation temperature, the sterilization time and culture media concentration on the time (days) to full colonization of the substrate for spawn multiplication by using the central composite design, the following procedure was adopted.

Three levels of PDA extract agar concentration (35 g, 50 g and 65 g) were weighed separately using laboratory analytical balance. Each concentration was dissolved in 1000 ml of distilled water and boiled for 2 minutes to dissolve the media. The media solution was then autoclaved for 10, 15 and 20 minutes each time but separately. This was followed by aseptically pouring of the cooled media into petri dishes in the laminar air flow which had been thoroughly sterilized using 70% alcohol. Once the media solidified in the Petri dishes, the oyster mushroom (Pleurotus ostreatus var. florida) was inoculated. This was achieved by cutting 2 × 2 mm² pieces of pure mycelia which was centrally placed on the cooled media. Each concentration was replicated 5 times making 15 plates for the three autoclaving and incubation temperatures per batch. The plates were completely sealed with a parafilm and incubated in dark condition at 20˚C, 30˚C and 25˚C each separately until the mycelium developed and covered the full area or otherwise. A complete randomized design was used to place the Petri plates in the incubator. The mycelia of Pleurotus ostreatus species were observed daily and measurements of the colony diameter of mycelia was noted after inoculation using a clear (transparent) ruler until the plates were fully colonized or otherwise.

The Petri dishes were regularly inspected and the contaminated ones were immediately discarded. The inset (b) of Figure 2 shows a contaminated Petri dish while the inset (d) is a fully colonized Petri dish. The mycelia growth covered the entire substrate within the 10 - 20 days of inoculation, by which the spawns were ready for use.

4.1.2. Analysis of a First-Order Model

Method of steepest descent was used to establish a point of local minimum, where the gradient was converging to zero. The level of temperature, time of sterilization and the culture media concentration level were sequentially altered along the path of steepest descent until there was no further decrease of days to full colonization.

Given that the Response Surface Methodology (RSM) is a sequential procedure

the first-order model $\stackrel{⌢}{y}={\stackrel{⌢}{\beta }}_0+{\displaystyle \underset{i=1}{\overset{k}{\sum }}{\stackrel{⌢}{\beta }}_i}{x}_i$ was fitted to lead rapidly and efficiently

along the path of improvement towards the general vicinity of the optimum. The curvature occurrence prompted the application of a second order model.

4.1.3. Analysis of a Second-Order Model

A model that incorporated the curvature was appropriate to approximate the response, given that the region of optimization was identified, the second order model was fitted to determine the optimum set of operating levels for the temperature, time of sterilization and the culture media concentration that minimized the days to full mycelium colonization.

4.1.4. Rotatable Second-Order Designs

To make a design both (approximately) orthogonal and rotatable, the axial distance for rotatability should be chosen, and then the center points are added [18]. Rotatability is a desirable and necessary property in response surface methodology for second order models exploration [16].

Under this design criterion, $Var\left[\frac{\partial \hat{y}\left(x\right)}{\partial x_i}\right]$ was constant $\forall x\in X$ that were

equidistant from the design center ( $i=1,2,\cdots ,k$ ). Therefore by applying Equation (1) the experimental runs were 19 and it was possible to have a rotatable and orthogonal design by adding the star- points and centre points replicated five times to the simple square factorial design points;

Such that given k = 3, p = 0, $n_f=n_s=1$ and $n_c=5$ using Equation (5) for the design to be rotatable it had axial distance (alpha) determined as;

$\alpha ={\left(\frac{2^{3-0}\left(1\right)}{1}\right)}^{\frac{1}{4}}\iff 8^{\frac{1}{4}}=\pm 1.682$

Since the value of α for a rotatable CCD depended on the fourth root of the cube [19].

To simplify the calculations and for the convenience sake the region of interest was defined by simple lower and upper limits on each of the design variables and coded as (−1, 1) interval hence making all the variables to be bound in the cube, $-1\le x_i\le 1$, with the centre points being the midpoint of these two values and assigned to “0”.

By letting ${\xi }_1$, ${\xi }_2$ and ${\xi }_3$ be the natural variables of temperature level, sterilization time and culture media concentration levels respectively, then the provided datasets used were factor settings of $\text{Temperature}=25\pm 5$,

$\text{SterilizationTime}=15\pm 5$ and the $\text{Culture}\text{media}\text{concentration}\text{level}=50\pm 15$ with five center point. Thus the coded variables for the temperature level, ster-

ilization time and culture media concentration were derived as; $x_1=\frac{{\xi }_1-25}{5}$, $x_2=\frac{{\xi }_2-15}{5}$ and $x_3=\frac{{\xi }_3-50}{15}$ respectively for temperature, as per Equation (7).

The actual data value for each input variable was sequentially altered and noted, as the process and the response surface was observed. The length of the colony diameter of mycelia was noted daily and the time to full colonization was observed as summarized in Table 1.

The rows in Table 1 represented the runs while the columns represented the levels of the factors. The first row was run at the “low” level of all of the factors, the second row was run at the “high” level of temperature as a factor but the “low” levels of sterilization time and agar media concentration.

Although the variable names remained $x_1$, $x_2$ and $x_3$ and design order 1 - 19, the actual experiment implementation was randomized. All the necessary properties for CCD second-order (quadratic) model such as rotatability and orthogonality were observed. To make an unbiased estimate of pure error, the CCD comprised of five centre points.

Alteration of the temperature level, sterilization time and the culture media concentration that shortened the mycelia colonization time (days), were as summarized in Table 2.

The shortest time observed was six (6) days with nineteen (19) days being the longest time. By using the coded second model data, it was possible to get the design matrix (X), the model equation and carry out the computation of the criteria and their efficiencies as discussed in the following section.

4.1.5. Optimality Criteria

The general three factors rotatable CCD was modified by varying the number of

times the center point was replicated. The optimality criteria D-, A-, and E which are best suited for the model parameters estimation were chosen to test the effect of varying the number of centre points in the design.

Given the information matrix ( $X^{\prime }X$ ), it was possible to determine the design optimality criterion at varied number of center points.

$X^{\prime }X=\left(\begin{array}{ccccccccc}13.66& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 13.66& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 13.66& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 8& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 8& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 8& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 24.01& 8& 8\\ 0& 0& 0& 0& 0& 0& 8& 24.01& 8\\ 0& 0& 0& 0& 0& 0& 8& 8& 24.01\end{array}\right)$

By using A-criterion, whereby a design is said to be A-optimal if it minimizes the sum of the main diagonal element (trace) of ${\left(X^{\prime }X\right)}^{-1}$ denoted as $\text{tr}\left({\left(X^{\prime }X\right)}^{-1}\right)$. The sum of the diagonal was 0.7445928 in all the designs. The value changed only due to the number of experimental runs.

With regard to the D-optimal, the value that maximized $\left|X^{\prime }X\right|$ was 7.48 and it was constant among different number of design centre points. Finally, regarding the E-optimal the value that minimized the maximum eigenvalue of ${\left(X^{\prime }X\right)}^{-1}$, denoted as $\min {\lambda }_{\max }{\left(X^{\prime }X\right)}^{-1}$, was 0.125, and it was constant among the three different chosen center points. The findings for different criteria are summarized in Table 3.

Therefore the addition of center points on CCD second-order designs had no significant effects on the optimality. However, the five center points were chosen to improve in making an unbiased estimate of pure error, since for an unbiased estimate the CCD should comprise of three to five centre points [22].

4.1.6. The Second Order Model Regression Equation

The second degree order model was computed using the R commands and the outputs were as given in Tables 4-7.

Therefore based on the second order model regression output, the model Equation (8) with the fitted regression coefficients becomes Equation (15)

$\begin{array}{c}\hat{y}=7.16-1.20x_1-1.37x_2-1.86x_3+2.08x_1^2+1.20x_2^2+1.20x_3^2\\ +0.25x_1{x}_2+0.00x_1{x}_3+0.25x_2{x}_3\end{array}$ (15)

Significant Predictors

Based on the t and the p-values it was observed that the entire predictor variables were statistically significant for influencing the time to spawn full development but at a varied levels of significance. Temperature was only significant at 5% but sterilization time and media culture concentration at 1% significance level. The temperature quadratic term was statistically significant at 1% significance

Signif. codes: 0 “***” 0.001 “**” 0.01 “*” 0.05 “.” 0.1 “ ” 1; Multiple R-squared: 0.8075, Adjusted R-squared: 0.615; F-statistic: 4.195 on 9 and 9 DF, p-value: 0.02201.

level. The coefficient of determination R² indicated that the model could approximate the data on the design points at 80.75% (after adjusting 61.5%).

Table 5 outputs display the first-order models (FO), which specifies the first-order response surface (i.e., a linear function), the canonical analysis of the

response surface with two-way interaction model (TWI) and the pure quadratic terms. The table also included lack-of-fit test and the pure error values.

Since the lack of fit was statistically insignificant ( $\mathrm{Pr}=0.078>0.05$ ) then the model could predict the response variable appropriately.

The response surface stationary points output was in coded form, which could be transformed into their natural value through the equations;

$x_1=\frac{{\xi }_1-25}{5}$, $x_2=\frac{{\xi }_2-15}{5}$ and $x_3=\frac{{\xi }_3-50}{15}$

such that;

${\xi }_1=26.29˚\text{C}$, ${\xi }_2=17.36\text{minutes}$ and ${\xi }_3=60.95\text{g}/\text{L}$

All the eigenvalues were positive implying the stationary points were minimum optimising regions.

The information in Tables 4-7 was computer generated. However, the same results could be arrived at through manual calculations as discussed in the following section.

4.1.7. Canonical Analysis

The fitted quadratic model provided a noticeable response variation with clarity and the canonical analysis gave the coordinates of the estimated stationary points and the canonical directions from the points.

The fitted canonical form model is characterized by $\hat{y}\left(w_i\right)={\hat{y}}_s+{\lambda }_1{w}_1^2+{\lambda }_2{w}_2^2+\cdots +{\lambda }_k{w}_k^2$ Where the fitted value at the stationary point was ${\hat{y}}_s$, then ${\lambda }_j$ were the eigenvalues obtained by getting the roots of the determinant equation $\left|B-\lambda I\right|=0$ and $w_j^2\text{,}\left(j=1,\cdots ,k\right)$ were the eigenvectors.

With reference to the second degree model which can be transformed into matrix form as;

$X=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]$, $b=\left[\begin{array}{c}-1.20\\ -1.37\\ -1.86\end{array}\right]$ and $\hat{B}=\left[\begin{array}{ccc}2.08& 0.125& 0.00\\ 0.125& 1.20& 0.125\\ 0.00& 0.125& 1.20\end{array}\right]$

with $b_0=7.16$

The eigenvalues and the eigenvectors output calculated using the R commands were as follows;

${\lambda }_1=2.10$, ${\lambda }_2=1.32$ and ${\lambda }_3=1.07$

$V_1=\left[\begin{array}{c}0.990\\ 0.141\\ 0.020\end{array}\right]$, $V_2=\left[\begin{array}{c}-0.110\\ 0.694\\ 0.731\end{array}\right]$ and $V_3=\left[\begin{array}{c}-0.090\\ 0.726\\ -0.682\end{array}\right]$

Now note the following operations using the transformed matrix values, eigenvalues and the eigenvectors;

$\hat{B}{V}_1={\lambda }_1{V}_1=\left[\begin{array}{c}2.077\\ 0.295\\ 0.042\end{array}\right]$, $\hat{B}{V}_2={\lambda }_2{V}_2=\left[\begin{array}{c}-0.145\\ 0.886\\ 0.961\end{array}\right]$ and $\hat{B}{V}_3={\lambda }_3{V}_3=\left[\begin{array}{c}-0.963\\ 0.777\\ 0.730\end{array}\right]$

Now suppose we let

$p=\left[\begin{array}{c}{V^{\prime }}_1\\ {V^{\prime }}_2\\ {V^{\prime }}_3\end{array}\right]=\left[\begin{array}{ccc}0.990& 0.141& 0.020\\ -0.110& 0.674& 0.731\\ -0.090& 0.726& -0.682\end{array}\right]$

so that,

$p^{-1}=\left[\begin{array}{ccc}0.989& 0.141& 0.019\\ -0.111& 0.673& 0.731\\ -0.089& 0.725& -0.682\end{array}\right]$

then

$p{p}^{\prime }=p{p}^{-1}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$

Therefore p is an orthogonal matrix since its columns and rows are orthogonal unit vectors (orthogonal vectors).

Suppose we let

$\Lambda =\left[\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right]=\left[\begin{array}{ccc}2.10& 0& 0\\ 0& 1.32& 0\\ 0& 0& 1.07\end{array}\right]$

Such that

$x=Pw\Rightarrow \left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{ccc}0.990& -0.110& 0.090\\ 0.141& 0.674& 0.726\\ 0.020& 0.731& -0.682\end{array}\right]\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right]$

implying

$\begin{array}{l}w=P^{-1}x\\ \Rightarrow \left[\begin{array}{c}{w}_1\\ w_2\\ w_3\end{array}\right]=\left[\begin{array}{ccc}0.990& -0.110& -0.090\\ 0.141& 0.673& 0.725\\ 0.019& 0.731& -0.682\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}0.990x_1-0.110x_x-0.090x_3\\ 0.141x_1+0.673x_2+0.725x_3\\ 0.019x_1+0.731x_2-0.682x_3\end{array}\right]\end{array}$

Now by using the preceding equation definitions of x and w, Equation (16) was obtained

$\begin{array}{l}\text{Giventhat}{w}^{\prime }\Lambda w={\lambda }_1{w}_1^2+{\lambda }_2{w}_2^2+{\lambda }_3{w}_3^2\\ =\left[\begin{array}{ccc}0.990x_1+0.141x_2+0.020x_3& -0.110x_1+0.674x_2+0.731x_3& -0.090x_1+0.726x_2-0.682x_3\end{array}\right]\\ \left[\begin{array}{ccc}2.10& 0& 0\\ 0& 1.32& 0\\ 0& 0& 1.07\end{array}\right]\left[\begin{array}{c}0.990x_1+0.141x_2+0.020x_3\\ -0.110x_1+0.674x_2+0.731x_3\\ -0.090x_1+0.726x_2-0.682x_3\end{array}\right]\\ =2.10{\left(0.990x_1+0.141x_2+0.020x_3\right)}^2+1.32{\left(-0.110x_1+0.674x_2+0.731x_3\right)}^2\\ +1.07{\left(-0.090x_1+0.726x_2-0.682x_3\right)}^2\\ =2.08x_1^2+1.20x_2^2+1.20x_3^2+0.25x_1{x}_2+0.00x_1{x}_3+0.25x_2{x}_3\end{array}$ (16)

Implying therefore

$w^{\prime }\Lambda w=2.08x_1^2+1.20x_2^2+1.20x_3^2+0.25x_1{x}_2+0.00x_1{x}_3+0.25x_2{x}_3\simeq x^{\prime }\hat{B}x$

Therefore $\hat{y}\simeq x^{\prime }\hat{B}x\simeq w^{\prime }\Lambda w$ and in canonical form $\hat{y}={\lambda }_1{w}_1^2+{\lambda }_2{w}_2^2+{\lambda }_3{w}_3^3$

Then given that all the eigenvalues of B were positive, then the quadratic form $x^{\prime }\hat{B}x$ was positive definite, hence the stationary point of the response surface was a point of minimum.

4.1.8. Stationary Points

To determine the levels of temperature ( $x_1$ ), sterilization time ( $x_2$ ) and culture media concentration level ( $x_3$ ) that minimized time to full colonization of the agar medium, the partial derivatives defined in Equation (10) were conducted such that;

$\partial \stackrel{⌢}{y}/\partial x_1=\partial \stackrel{⌢}{y}/\partial x_2=\partial \stackrel{⌢}{y}/\partial x_3=0$

That is the stationary point were obtained by $\frac{\text{d}\hat{y}}{\text{d}x}=b+2\hat{B}x=0$

$\Rightarrow x_s=-\frac{1}{2}{B}^{-1}b$

and

${\hat{y}}_s=b_0+\frac{1}{2}{x^{\prime }}_{s}b$

Therefore given

$\hat{B}=\left[\begin{array}{ccc}2.08& 0.125& 0.00\\ 0.125& 1.20& 0.125\\ 0.00& 0.125& 1.20\end{array}\right],{\hat{B}}^{-1}=\left[\begin{array}{ccc}0.484& -0.051& -0.005\\ -0.051& 0.848& -0.088\\ 0.005& -0.088& 0.843\end{array}\right]$

and

$b=\left[\begin{array}{c}-1.20\\ -1.37\\ -1.86\end{array}\right]$

Then

$X_s=-\frac{1}{2}\left[\begin{array}{ccc}0.484& -0.051& -0.005\\ -0.051& 0.848& -0.088\\ -0.005& -0.088& 0.843\end{array}\right]\left[\begin{array}{c}-1.20\\ -1.37\\ -1.86\end{array}\right]=\left[\begin{array}{c}0.260\\ 0.480\\ 0.726\end{array}\right]$

and therefore

$x_1=\frac{{\xi }_1-25}{5}\Rightarrow 0.260=\frac{{\xi }_1-25}{5}$ $\therefore {\xi }_1=26.3$

$x_2=\frac{{\xi }_2-5}{5}\Rightarrow 0.480=\frac{{\xi }_2-15}{5}$ $\therefore {\xi }_2=17.4$

$x_3=\frac{{\xi }_3-50}{15}\Rightarrow 0.726=\frac{{\xi }_3-50}{15}$ $\therefore {\xi }_3=60.89$

The computer generated stationary points output figures and the manually computed were exactly the same when rounded to one decimal place as shown in the summary findings of the manually computed values in Table 8.

Therefore the response surface could be predicted by using the following equation

${\hat{y}}_s=b_0+\frac{1}{2}{x^{\prime }}_{s}b$

such that

${\hat{y}}_s=7.16+\left[\begin{array}{ccc}0.260& 0.480& 0.726\end{array}\right]\left[\begin{array}{c}-1.20\\ -1.37\\ -1.86\end{array}\right]=7.16-1.152=6.008\text{Days}$

4.1.9. Diagnostics and Plots of the Estimated Response Surfaces

The basic assumptions on errors which include, independence, normality and the constancy of variance for errors in order to apply the second order statistical model and ANOVA test statistics were tested. The findings are summarized in Figure 3.

The fact that there were equally spread residuals around the horizontal line without any distinct patterns in the residual versus fitted values plot, then that was an indication of residuals linear or non-linear relationships nonexistence, thereby satisfying the assumption of residuals independence. The normal quantile-quantile plot shows the residuals were normally distributed given that the observed values were along the reference line hence ascertained normality. The scale-Location plot showed that the residuals were equally spread along the ranges of predictors and hence the assumption of equal variance (homoscedasticity) was not violated. Based on the residuals versus leverage plot, most of the predictor variables were well inside of the Cook’s distance lines hence few identified outliers. The following response surface graphs display the predictor variables

on the x- and y-axes and then a continuous surface that represents the fitted response values on the z-axis. For each of the surface map, two variables were displayed at a time holding the third variable constant.

Figure 4 shows the response surface map and a contour plot for temperature and time of substrate sterilization for spawn’s multiplication.

They represent a minimum response surface thus indicating a point of optimum operating conditions reaching a minimum. By inspection the plots are non-linear, implying a strong temperature and sterilization interaction effect on the time to spawns maturity and the response surface is curved due to the fact that the model contains quadratic terms that were statistically significant.

The contour lines appeared to follow the direction of movement along the path of minimum response from a reference point. The contour lines were within the response values in the data which ranged between 6 to19 nevertheless it is worth noting that the accuracy of the response surface plot depends on how well the model represents the true relationships among the variables.

The lowest values of days to full colonization of the substrate in the petri dish were in the upper side of the contour plot as it appears in Figure 5, which corresponds with high values of both agar media concentration and temperature, time of sterilization was held constant for this plot.

As with the contour plot the lowest values of the response surface were along the side opposite the temperature axis, which corresponded with high values agar media concentration and temperature

The following and the last is the contour plot and then the response surface map for agar media concentration and time of sterilization.

As in the contour plot in Figure 6, the lowest values of days to full substrate colonization are in the lower right corner of the response surface, which corresponded with high values of agar media concentration.

The lowest values of days to full substrate colonization are in the upper right corner of the plot, which corresponds with high values of agar media concentration but average time of sterilization, temperature was held constant in this particular analysis and representation.

5. Conclusion and Recommendations

The CCD provided the statistical elements necessary for the evaluation of Temperature, agar concentration and sterilization time during spawn propagation. The contour plots and the response surface maps generated characterized surface and helped a great deal in locating the optimum points precisely at a glance. Routinely training on how one would process the spawns using the locally available agricultural by-products or forest wastes, using the scientifically proven formulations was timely and would not only provide a protein-rich food but reduces the environmental pollution. It is a transformative link for inedible wastes into edible biomass of high monetary value.

Multiple response optimizations aim to achieve maximal nutritional value and yield against minimal duration of spawns propagation and the spawns quality is suggested and recommended.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Anderson, C.M., Montgomery, C.D. and Myers, H.R. (2016) Response Surfaces Methodology: Process and Product Optimization Using Designed Experiments. 4th Edition, Wiley, Hoboken, 111, 297-320.
[2]	Diamond, W.J. (1981) Practical Experiment Design for Engineers Scientists. Lifetime Learning Publications, New Jersey, USA.
[3]	Chang, S.T. (1996) Mushroom Research and Development-Equality and Mutual Benefit. In: Royse, D.J. (Ed.), Mushroom Biology and Mushroom Products, Proceedings of the 2nd International Conference, Pennsylvania State University, USA, 1 September 1996, 473-476.
[4]	Wambua, J. (2014) Mushroom Cultivation in Kenya. In: Mushroom Growers Handbook, Oxford University Press, Nairobi, 197-123.
[5]	Sharma, V.P., Kumar, R., Gupta, R.K., Satish, K. and Rajender, S. (2013) Optimizations of Parameters for Quality Spawn Production. Mushroom Research, 21, 31-36.
[6]	Ansari, A.N. and Elhami, B. (2008) Effect of Substrates of Spawn Production on Mycelium Growth of Oyster Mushroom Species. Journal of Biological Sciences, 8, 474-477. https://doi.org/10.3923/jbs.2008.474.477
[7]	Samuel, A.A. and Eugene, T.L. (2012) Growth Performance and Yield of Oyster Mushroom (Pleurotus ostreatus) on Different Substrates Composition in Buea South West Cameroon. Science Journal of Biochemistry, 2012, sjbch-139. https://doi.org/10.7237/sjbch/139
[8]	Hoa, H.T. and Wang, C.-L. (2015) The Effects of Temperature and Nutritional Conditions on Mycelium Growth of Two Oyster Mushrooms (Pleurotus ostreatus and Pleurotus cystidiosus). Mycobiology, 43, 14-23. https://doi.org/10.5941/MYCO.2015.43.1.14
[9]	Rajmohan, T. and Palanikumar, K. (2012) Application of the Central Composite Design in Optimization of Machining Parameters in Drilling Hybrid Metal Matrix Composites. Measurement, 46, 1470-1481. https://doi.org/10.1016/j.measurement.2012.11.034
[10]	Al-Shingiti, A.M. and Huda, S. (2004) On Second-Order A-, D- and E-Minimax Designs for Estimating Slopes in Extrapolation and Restricted Interpolation Regions. Communications in Statistics—Simulation and Computation, 33, 773-785. https://doi.org/10.1081/SAC-200033245
[11]	Bahrim, G., Horincar, G., Popa, A. and Vincetiu (2017) Submerged Cultivation of Fomes fomentarius Mushroom and Increase of Biomass Yield by Statistical Design of Experiments and Mathematical Modeling. Dunarea de Jos University of Galati, Faculty of Food Science and Engineering, Galati.
[12]	Dressaire, E. (2016) Mushrooms Use Convectively Created Airflows to Disperse Their Spores. Proceedings of the National Academy of Sciences of the United States of America, 113, 2833-2838. https://doi.org/10.1073/pnas
[13]	Sarrai, A.E., Hanini, S., Merzouk, N.K., Tassalit, D., Tibor, S., Hernádi, K. and László Nagy, L. (2016) Using Central Composite Experimental Design to Optimize the Degradation of Tylosin from Aqueous Solution by Photo-Fenton Reaction. Materials, 9, 428. https://doi.org/10.3390/ma9060428
[14]	Li, P.Q., et al. (2012) Medium Optimization for Exopolysaccharide Production in Liquid Culture of Endophytic Fungus Berkleasmium sp. Dzf12. International Journal of Molecular Sciences, 13, 11411-11426. https://doi.org/10.3390/ijms130911411
[15]	Ahmadi, M., Vahabzadeh, F., Bonakdarpour, B., Mofarrah, E. and Mehranian, M. (2005) Application of the Central Composite Design and Response Surface Methodology to the Advanced Treatment of Olive Oil Processing Wastewater Using Fenton’s Peroxidation. Journal of Hazardous Materials, 123, 187-195. https://doi.org/10.1016/j.jhazmat.2005.03.042
[16]	Box, G.E.P. and Hunter, J.S. (1957) Multi-Factorial Experimental Designs for Exploring Response Surfaces. Annals of Mathematical Statistics, 28, 195-241. https://projecteuclid.org/euclid.aoms/1177707047 https://doi.org/10.1214/aoms/1177707047
[17]	Box, G.E.P. and Draper, N.R. (1987) Wiley Series in Probability and Mathematical Statistics. Empirical Model-Building and Response Surfaces. John Wiley & Sons, Hoboken.
[18]	Cornell, J.A. and Khuri, A.L. (1987) Response Surfaces: Designs and Analyses. Marcel Dekker Inc., New York, 405.
[19]	Montgomery, D.C. (1991) Design and Analysis of Experiments. John Wiley and Sons, New York, 204.
[20]	Draper, N.R. (1992) Introduction to Box and Wilson (1951) on the Experimental Attainment of Optimum Conditions. In: Kotz, S. and Johnson, N.L., Eds., Breakthroughs in Statistics, Springer Series in Statistics (Perspectives in Statistics), Springer, New York, 267-269. https://doi.org/10.1007/978-1-4612-4380-9_22
[21]	Lenth, R.V. (2012) Response-Surface Methods in R, Using rsm. The University of Iowa, Iowa City, 7-11. https://cran.r-project.org/web/packages/rsm/vignettes/rsm.pdf
[22]	Chigbu, P.E. and Orisakwe, N.C. (2000) On Optimal Partially Replicated Rotatable and Slope Rotatable Central Composite Designs. Department of Statistics University of Nigeria, Nsukka.