Introduction of a mathematical model for optimizing the drug release in the patient’s body

Abstract

Background

Drug release in a patient’s body is of particular interest to the pharmaceutical industry. One of the most essential types of drug release is the gradual release based on a behavior, which is called a profile or modified release. The investigation of the time-oriented quality characteristic is one of the newest topics in the area of product design. There are already several approaches addressing this issue. In this paper, a mathematical model is proposed to find the suitable values of the controllable factors in a drug to achieve the profile of the drug release in the patient’s body.

Results

The proposed method has several advantages over the existing methods.

Conclusion

The authors feel that by adjusting the control factors during the production process the drug release profile become closer to the reference profile.

Introduction

The amount of time it takes a drug to release in a patient’s body as well as the time it takes to exert its effects on the target organ are very important factors used to measure the effectiveness of a drug. If this releasing manner is not based on a pre-defined profile, it may cause a reduction of curative properties of the drug and can even have some negative effects on the patient’s body. Similarly, in the area of quality engineering, the time-oriented quality characteristics are also assessed. The time-oriented profile of the quality characteristic is specified and the aim of the designer is to find the predefined profile with minimum deviation from the target. The quality characteristics are then monitored using the defined profile. In this study, we aim to establish a logical relationship between these two areas and to apply a mathematical modeling approach to investigate the drug release problem in pharmaceutics. In this paper some basic definitions of drug release and quality engineering are presented and then we introduce the four existing approaches for these types of problems and their deficiencies. The proposed method is presented in the next section. Several examples are provided to evaluate the suggested model and in the final section, the conclusions are made.

Definitions

In this section some of the basic terms included in the paper are defined to familiarize the reader with the concepts of the discussion.

Drug release

Drug release is an important stage in the drug life cycle. When the drug is released based on a pre-defined profile, it is more effective on the patient’s body. One of the most applicable approaches for measuring the amount of released drugs is to measure the plasma concentration of the drug. The drug is considered effective when the plasma concentration is somewhere between minimum effective concentration (MEC) and minimal toxic concentration (MTC) as is shown in Figure 1[1–3].

Figure 1

Plasma concentration versus time profile[1].

Drugs are usually classified based on the drug release mechanism as follows:

Immediate release drugs: In this group, the drug is quickly released in the body. This is particularly suitable for drugs that need to take affect rapidly such as painkillers [1, 4].

Modified drug release: In this case by using the pharmaceutical techniques, the time, the amount and the target organ for the drug release is determined. The delayed release and extended release are the methods being used. In the delayed release the drug is released after a pre-determined delay. Figure 2 shows the plasma concentration for this modified release method [1, 4].

Figure 2

Plasma concentration versus time profile for an immediate release drug and a delayed release drug[1].

In the extended release technique, the drug is released gradually over a longer period of time. It is classified into two categories: sustained release and controlled release. In sustained release, the drug is released continuously with a constant rate. In controlled release, the drug is released intelligently so that the concentration remains almost constant in the body. Figure 3 shows the plasma concentration when using this method of drug release [1, 5].

Figure 3

Plasma concentration versus time profile for a controlled release formulation[1].

Time-oriented quality characteristics

There are several definitions of the quality characteristics in the quality management literature. The most comprehensive of them is the degree of adaptability of the quality characteristic by the user’s requirements [6]. Furthermore, the design phase is the principal stage of a product life cycle, because the quality is formed in this stage and control actions at the end of the production process cannot improve the quality of a product with poor quality of design [7].

The Taguchi robust design is a famous design procedure. It is an engineering method for optimizing the product or process condition to minimize the product sensitivity to the noise factors in the environment, such as: ambient temperature, humidity, air pressure and direct sunlight [8]. So, a product with high quality and low cost is being produced. One property of this approach is to investigate the quality characteristics numerically. In this approach the quality characteristics are grouped into three classes as: nominal the best (NTB), larger the better (LTB) and smaller the better (STB). Each of these quality characteristics could be constant or variable over time [9].

The target value and the specification limits for the time-oriented quality characteristics are being changed over time. So, for the design of a product with these quality characteristics, the parameters are designed such that the quality characteristics are being as close to their pre-specified target values as possible.

In this regard, three basic topics need to be introduced.

Design of experiments (DOE)

A collection of statistical methods that are used to find the influenced factors on a quality characteristic and to optimize its conditions. There are several types of DOE techniques including factorial experiments and fractional factorial experiments [10, 11].

Response surface methodology (RSM)

A statistical and mathematical method for modeling, analyzing and optimizing the problems with response variables which are directly related to some other independent variables [12].

Desirability function

Is one of the common methods to simultaneously optimize multi response problems. The most applicable method of this type is the Derringer and Suich’s which is defined for several types of quality characteristics as follows [13]:

NTB quality characteristic:

$$\mathit{DF}\left(y\right)=\left\{\begin{array}{c}\hfill {\left[\frac{y-\mathit{LSL}}{T-\mathit{LSL}}\right]}^r,\mathit{LSL}<y<T\hfill \\ \hfill {\left[\frac{y-\mathit{USL}}{T-\mathit{USL}}\right]}^s,T<y<\mathit{USL}\hfill \\ \hfill 0,y<\mathit{LSL};y>\mathit{USL}\hfill \end{array}\right.$$

(1)

LTB quality characteristic:

$$\mathit{DF}\left(y\right)=\left\{\begin{array}{c}\hfill 1,y>y_i^{*}\hfill \\ \hfill {\left[\frac{y-y_{i_{*}}}{y_i^{*}-y_{i_{*}}}\right]}^r,y_{i_{*}}<y<y_i^{*}\hfill \\ \hfill 0,y<y_{i_{*}}\hfill \end{array}\right.$$

(2)

STB quality characteristic

$$\mathit{DF}\left(y\right)=\left\{\begin{array}{c}\hfill 1,y<y_{i_{*}}\hfill \\ \hfill {\left[\frac{y_i^{*}-y}{y_i^{*}-y_{i_{*}}}\right]}^r,y_{i_{*}}<y<y_i^{*}\hfill \\ \hfill 0,y>y_i^{*}\hfill \end{array}\right.$$

(3)

In the above equations:

y: value observed for the quality characteristic

T: The target value for quality characteristic applicable for NTB quality characteristic.

USL: Upper specification Limit of NTB quality characteristic

LSL: Lower specification Limit of NTB quality characteristic

y _i *: optimum point for LTB quality characteristic and highest acceptable value for STB quality characteristic

y _i* : Optimum point for STB quality characteristic and lowest acceptable value for LTB quality characteristic

r, s: Weight values, positive constants.

Problem definition

The drugs have a pre-determined profile for release based on the drug’s controlled-release mechanism. The aim in any drug laboratory is to find optimum adjustment of the controllable factors, such as material, production machine settings and so on to produce drugs that achieve the predetermined profile as much as possible. Four methods already exist for parameter design of a drug to achieve its pre-determined profile:

Contour overlay method

This method is applied by Gohel and Amin [14] to find the optimal values to the Diclofenac Sodium formulation. The aim is to determine the suitable values for the three main controllable factors: stirring speed, concentration of CaCl₂ and percentage of liquid paraffin, all of which influence the drug efficacy. The pre-determined profile of release is defined in advance. Then, the regression function of the drug release as a response variable and the above-mentioned control factors as independent variables is obtained by the least square method. For each point of time, the response is computed and compared to the pre-specified value. In this method, one variable is kept fixed and a two dimensional plot is used to find the optimal values.

The disadvantage of this method is that when the number of control factors increases, the efficiency of the method to introduce optimal values decreases.

Profile selection

In situations where the profile properties are hard to identify, selection of the best profile is done by using the pre-defined indices. Two of these indices are f₁ and f₂ defined as:

$$f_1=\frac{{\displaystyle {\sum }_{t=1}^n}\left|R_t-T_t\right|}{{\displaystyle {\sum }_{t=1}^n}{R}_t}*100$$

(4)

$$f_2=50*\mathit{Log}\left\{{\left[1+\frac{1}{n}{\displaystyle \sum _{t=1}^n}{\left(R_t-T_t\right)}^2\right]}^{-0.5}*100\right\}$$

(5)

Where:

R_t: Percentage of drug release obtained from the reference formulation

T_t: Percentage of drug release obtained from the test formulation

n: number of observations

The first index, f₁ is defined as the dissimilarity index. As long as its value is small; the profile is close to the reference profile. The second index, f₂ is defined as similarity index and when its value is large; the profile is near to the reference profile [15, 16].

MSE minimizing method

This method is applied in three articles. Truong et al. [17] used this method to determine the optimum values for control factors of a regenerative drug based on a profile of seven points.

Park et al. [18] used this method to investigate two quality characteristics separately for six and seven point profiles. Shin et al. [19] used this method to assess two quality characteristics separately for eight and eleven point profiles.

The first step in this method is to gather data and to calculate the basic statistics such as the mean and the variance. Then the RSM for these statistics are computed at each point of time. The optimal values for the control factors are obtained such that the following objective function is minimized.

$$\mathit{Minimize}{\displaystyle \sum _{q=1}^w}{\left(\widehat{M}\left(x,t_q\right)-T_q\right)}^{\mathit{2}}+{\displaystyle \sum _{q=1}^w}\widehat{v}\left(x,t_q\right),S\mathit{.}t:x\in \Omega$$

(6)

Where:$\widehat{M}\left(x,t_q\right)$: The mean of the responses at time t_q.

$\widehat{v}\left(x,t_q\right)$: The variance of the responses at time t_q.

T_q: The pre-specified target value for the response variable for the time q.

w: The number of points in time under study.

Method of minimizing the total cost

This method is used by Goethals and Cho [20] and also the experiment of Gohel and Amin [14] on the Diclufenac Sodium is reassessed. The logic behind this method is to find the optimal values for control factors that minimize the following objective function:

$$\begin{array}{ll}\mathit{Minimize}E\left[\mathit{TC}\right]& ={\displaystyle \sum _{q=1}^w}[{\displaystyle \underset{\mathit{LS}{L}_q}{\overset{\mathit{US}{L}_q}{\int }}}L\left[Y\left(q\right)\right]\mathit{.}f[Y\left(q\right)]\mathit{.}\mathit{dY}\left(q\right)\\ +{\displaystyle \underset{-\infty }{\overset{\mathit{LS}{L}_q}{\int }}}N{C}_{q1}\mathit{.}f\left[y\left(q\right)\right]\mathit{dY}\left(q\right)\\ +{\displaystyle \underset{\mathit{US}{L}_q}{\overset{+\infty }{\int }}}N{C}_{q2}\mathit{.}f\left[Y\left(q\right)\right]\mathit{dY}\left(q\right)]\end{array}$$

(7)

Where:

LSL_q and USL_q: are the lower and the upper specification limits for the quality characteristic, respectively.

f(y(q)): is the probability distribution function for response variable at time q.

NC_q1 and NC_q2: are the costs corresponding to being greater than USL and smaller than LSL, respectively.

L(y (q)): is the quality loss function for the quality characteristic within the acceptable region, but not on the target.

w: is the number of time points under study.

The proposed method

The proposed method is a systematic and straightforward technique for determining the optimum values for the control factors for a drug. So that in a specified time interval, the drug release follows its premeditated profile. This method requires the following steps:

• 1. Determination of the drug release profile: Considering the kind of drug and its mechanism of release, the pharmaceutics design of the release profile of a drug by consulting the specialist physicians. To facilitate the comparison between the standard profile and the drug profile function, some points on time are considered and the experiments are run in these points. At each time point, the target value and the upper and the lower specification limits are determined. Selection of the number of points under study is based on the type of the drug and its life cycle in the patient’s body.

• 2. Determination of the experiment templates: In this stage, many controllable factors such as raw material and production factors for the drug under study are determined. Several combinations of these controllable factors are being tested by running the experiments. One important logic of the DOE is to find as much as information possible from the minimum number of experiments. For each combination of the factor levels at each time point some data is collected. Then, the data are organized based on the Table 1.The primary statistics such as the mean, the variance and the coefficient of variation for each time point and the covariance between observations in different time points are calculated. The computational formulas used to compute these statistics are as follows:

Table 1 Experimental format[20]

$${\overline{y}}_{\mathit{qr}}=\frac{{\displaystyle {\sum }_{w=1}^m}{y}_{\mathit{qrw}}}{m}$$

(8)

$$s_{\mathit{qr}}^2=\frac{{\displaystyle {\sum }_{w=1}^m}{\left(y_{\mathit{qrw}}-{\overline{y}}_{\mathit{qr}}\right)}^2}{m-1}$$

(9)

$${\left(\frac{s}{m}\right)}_{\mathit{qr}}=\frac{s_{\mathit{qr}}}{{\overline{y}}_{\mathrm{qr}}}$$

(10)

$$s_{i,j}=\frac{{\displaystyle {\sum }_{r=1}^m}\left(y_{\mathit{ipr}}-{\overline{y}}_{\mathit{ip}}\right)\left(y_{\mathit{jpr}}-{\overline{y}}_{\mathit{jp}}\right)}{m-1}$$

(11)

• 3. Determination of the relationships among the statistics and the control factors: By using RSM technique, the relationships are defined. For the sake of simplicity and prevention of using data with several scales, the control factors are coded by linear relationships.

  $$\begin{array}{l}{\widehat{\mu }}_q\left(x\right)=x{\widehat{\beta }}_{\mathit{\mu q}},{\widehat{\beta }}_{\mathit{\mu q}}={\left(x^{\text{'}}x\right)}^{-1}{x}^{\text{'}}{\overline{y}}_q,x=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \cdots \hfill & \hfill x_{1,k-1}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 1\hfill & \hfill \cdots \hfill & \hfill x_{n,k-1}\hfill \end{array}\right]\\ {\overline{y}}_q=\left[{\overline{y}}_{q1},{\overline{y}}_{q2},\dots ,{\overline{y}}_{\mathit{qn}}\right]\text{'}\end{array}$$

  (12)
  $$\begin{array}{l}{\widehat{s}}_q^2\left(x\right)=x{\widehat{\beta }}_{s^{2}q},{\widehat{\beta }}_{s^{2}q}={\left(x^{\text{'}}x\right)}^{-1}{x}^{\text{'}}{s}_q^2,x=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \cdots \hfill & \hfill x_{1,k-1}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 1\hfill & \hfill \cdots \hfill & \hfill x_{n,k-1}\hfill \end{array}\right]\\ s_q^2=\left[s_{q1}^2,s_{q2}^2,\dots ,s_{\mathit{qn}}^2\right]\text{'}\end{array}$$

  (13)
  $$\begin{array}{l}{\left(\left(\widehat{\frac{s}{m}}\right)\right)}_q\left(x\right)=x{\widehat{\beta }}_{\left(s/m\right)}q,{\widehat{\beta }}_{\left(s/m\right)}q=\left(x^{\text{'}}x\right){}^{\left(-1\right)}{x}^{\text{'}}(s/m){}_q,\\ x=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \cdots \hfill & \hfill x_{1,k-1}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 1\hfill & \hfill \cdots \hfill & \hfill x_{n,k-1}\hfill \end{array}\right]\\ {\left(\frac{s}{m}\right)}_q=\left[\right(\frac{s}{m}){}_{q1},(\frac{s}{m}){}_{q2},\dots ,(\frac{s}{m}\left){}_{\mathit{qn}}\right]{}^{\text{'}}\end{array}$$

  (14)
  $$\begin{array}{l}{\widehat{s}}_{i,j}\left(x\right)=x{\widehat{\beta }}_{s_{i,j}},{\widehat{\beta }}_{s_{i,j}}={\left(x^{\text{'}}x\right)}^{-1}{x}^{\text{'}}{s}_{i,j},x=\left[\begin{array}{ccc}\hfill 1\hfill & \hfill \cdots \hfill & \hfill x_{1,k-1}\hfill \\ \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \hfill 1\hfill & \hfill \cdots \hfill & \hfill x_{n,k-1}\hfill \end{array}\right]\\ s_{i.j}=\left[s_{i,j,1},s_{i,j,2},\dots ,s_{i,j,n}\right]\text{'}\end{array}$$

  (15)

In the interest of time and cost, the number of control factors is reduced before running the experiments by using any technique such as screening experiments, as well as the forward, backward and stepwise regression.

• 4 Model optimization: Using the desirability function method, the optimal values for control factors are determined based on the type of quality characteristics and their specification limits such that their values come as close to the target values as possible. The desirability function of interest is:

  $$\begin{array}{ll}\mathit{Maximize}{D}_{\mathit{total}}& =\left\{\left[{\displaystyle \prod _{i=1}^n}D{\left({\mu }_i\right)}^{w_i}\right].\left[{\displaystyle \prod _{i=1}^n}D{\left(s_i^2\right)}^{w_i^{\text{'}}}\right]\right.\\ \left..\left[{\displaystyle \prod _{i=1}^n}D{\left(\frac{s_i}{m_i}\right)}^{w_i^{\text{'}\text{'}}}\right].\left[{\displaystyle \prod _{i=1}^n}D{\left(s_{i,j}\right)}^{w_i^{\text{'}\text{'}\text{'}}}\right]\right\}\\ \times \left(\frac{1}{{\displaystyle \sum _{i=1}^n\left(w_i+w_i^{\text{'}}+w_i^{\text{'}\text{'}}\right)+{\displaystyle \sum _{i=1}^{\left(\begin{array}{l}n\\ 2\end{array}\right)}{w}_i^{\text{'}\text{'}\text{'}}}}}\right)\end{array}$$

  (16)

The results are robust as long as the covariances between the observations for each pair of points are close to zero. So, when there is a deviation in some time intervals, they would not be transmitted to the other points.

The other advantage of the proposed method is its ability to be used for any part of the desirability function. For instance when we don’t have access to the entire data and only the mean and the variance of the observations are available, the covariance part of the model may be eliminated. Or if the mean of the observations at each point of time for different combinations is in hand, only the mean part of the model is being used. Also, by using the desirability function and its weighted values, one may use any indices in some points under study. For the sake of simplicity, in the examples provided in Section 5, equal weights are assigned to all statistical indices in all time periods.

Numerical examples

To illustrate the applications of the proposed method, seven examples for different drugs are presented in this section adapted from credible pharmaceutical papers. These examples are solved by the proposed method to find the optimum values for the control factors of the drugs. The required material, the methods of pharmaceutical experiments and the data for each example are presented in the stated indicated references.

Example 1

Diclufenac Sodium

The release profile of this drug is investigated by Gohel [14] and Goethals [20]. The contour overlay and the minimization of quality loss function methods are introduced in their papers, respectively. This drug has three main control factors given in Table 2.

Table 2 Main control factors influencing Diclufenac Sodium release

The first step is to code the control factors using the following relationships:

$$x_{1\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_1-1000}{500},x_{2\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_2-10}{5},x_{3\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_3-25}{25}$$

In this research, three points of time for the drug release profile are being investigated with properties shown in Table 3.

Table 3 The target values and lower and upper values for example 1

The response surface relationships for the mean, the variance, the coefficient of variation and the covariance between each pair of points under study are presented in the Appendix 1. Optimum values are shown in Table 4.

Table 4 Optimum values for example 1

Example 2

Terazosin HCl dehydrate

The release profile for this drug is investigated by Shin [19] and the problem is solved by the MSE minimization method. This experiment has ten control factors as shown in Table 5.

Table 5 Control factors influencing Terazosin HCl dehydrate release

Noticing the large number of control factors in this example, five control factors x₁, x₃, x₇, x₈ and x₁₀ are identified as significant control factors by using the stepwise regression method. The control factors are coded by the following relationships:

$$x_{i-\mathit{new}}=\left\{\begin{array}{c}\hfill \frac{x_i-93.71}{7.03},i=1\hfill \\ \hfill \frac{x_i}{7.03},i=2,3,\dots ,10\hfill \end{array}\right.$$

In this research, 11 points of time of drug release profile are being investigated as presented in Table 6.

Table 6 The target values and lower and upper values for example 2

The response surface relationships for the mean and the variance of the underlying data are presented in Appendix 2. Optimum values for this example are shown in Table 7.

Table 7 Optimum values for control factors for example 2

Example 3

Verapamil HCl

The release profile of this drug is investigated by Siva [21]. The three main control factors for this drug are presented in Table 8.

Table 8 Main control factors influencing Verapamil HCl release

The control factors are coded by the following relationships:

$$x_{1\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_1-11}{3},x_{2\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_2-36}{12},x_{3\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_3-90}{30}$$

In this research, five points of time are investigated from release profile as shown in Table 9.

Table 9 The target values and lower and upper values for example 3

The RSM relationships for the mean, the variance and the coefficient of variation for the points in Table 8 are presented in Appendix 3. By using the desirability function method the optimum values obtained for control factors are shown in Table 10.

Table 10 Optimum values for control factors for example 3

Example 4

Metformin

The release profile for this drug is investigated by Nagrava [22]. The three main control factors are defined for this drug release as shown in Table 11.

Table 11 Main control factors for example 3

The values of the control factors are coded using the following relationships:

$$x_{1\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_1-1.758}{1.25},x_{2\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_2-0.25}{0.25},x_{3\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_3-3.75}{1.25}$$

The three points of time for the release profile are investigated in this research have the properties provided in Table 12.

Table 12 The target values and lower and upper values for example 4

The RSM relationships for the mean, the variance and the coefficient of variation for the data are presented in Appendix 4. By using the desirability function method the optimum values obtained for control factors are shown in Table 13.

Table 13 Optimum values of control factors for example 4

Example 5

Rhinetedin

The release profile of this drug is investigated by Patel [23]. The two main control factors for this drug are presented in Table 14.

Table 14 Main control factors for example 5

The control factors are coded by the following relationships:

$$x_{1\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_1-672}{168},x_{2\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_2-168}{84}$$

In this research three time points are investigated from release profile are shown in Table 15.

Table 15 The target values and lower and upper specifications for example 5

In this example, the index f₂ is the measure of similarity between the drug release profile and the target profile. The RSM relationships are presented in Appendix 5 and the optimum values are shown in Table 16.

Table 16 Optimum values for example 5

Example 6

Metoprolol

The release profile for this drug is investigated by Gohel [24]. The two main control factors defined for this drug are shown in Table 17.

Table 17 Main control factors for example 7

The control factor values are coded by using the following relationships:

$$x_{1\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_1-30}{10},x_{2\left(\mathrm{new}\right)}=\frac{{\mathrm{x}}_2-20}{10}$$

The three points of time for the drug release profile are presented in Table 18.

Table 18 The target values and lower and upper specification limits for example 7

In the study of this drug, f₂, t₅₀ (the time required for 50% of drug to be released) and mean dissolution time (MDT) are the measures of the similarity factor between release profile and the predefined profile, the time required to dissolve half of the drug and the mean dissolution time, respectively. The RSM relationships for the means and these measures are presented in Appendix 6. By using the proposed method, the optimum values are obtained as shown in Table 19.

Table 19 Optimum values for example 7

Comparison of the proposed method and the existing ones

The disadvantages of the existing methods are:

Contour overlay method:

This method has a limited application and when the number of variables exceeds from two, the model may not be optimized unless the additional variables are being fixed at a constant level.

Profile selection method:

In this method, the number of test profiles is adjusted based on the experimenter point of view and the best profile is selected among the existing ones. It is possible that the optimum values for the control factors may not be included in these profiles.

MSE minimizing method:

In this method, there is no attention paid to the specification limits, while in the real world, passing these limits has substantial penalties.

Minimizing the total cost method:

In this method all deviations from the target values are evaluated by means of money terms, while in human problems, e.g. pharmaceutical studies, adverse events may have human fallout which cannot be measured by money terms.

The proposed method overcomes all the above disadvantages.

Conclusions

Investigation of the pharmaceutics problems in an industrial engineering framework is very constructive. The key point here is the problem presentation by the engineering terms. In this research, the drug release problem which is an important subject of pharmaceutics is being studied. In this area, applying the complex formulas is avoided. So, the experts with minimum knowledge of mathematics and statistics may apply this approach to solve the pharmaceutics problems. The results of the examples show the ability of the proposed model for solving the controlled release problems and to assure that the intended drug is resolved as its predefined profile. The simultaneous optimization of drugs with multi time-oriented quality characteristics is a topic for the future research.

Appendix 1

$$\begin{array}{ll}{\mu }_{1\left(1h\right)}& =39.929+2.365x_1-2.206x_2-1.959x_3\\ +0.202x_1^2+1.971x_2^2-0.912x_3^2-1.389x_1{x}_2\\ +0.797x_1{x}_3+0.079x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{2\left(6h\right)}& =73.368+4.388x_1-5.031x_2-2.379x_3\\ +0.399x_1^2+0.579x_2^2-0.127x_3^2-1.525x_1{x}_2\\ -0.062x_1{x}_3-0.359x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{3\left(8h\right)}& =83.203+4.165x_1-4.562x_2-2.498x_3\\ -0.624x_1^2-0.907x_2^2+1.176x_3^2-2.37x_1{x}_2\\ +0.151x_1{x}_3-1.632x_2{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{1\left(1h\right)}& =7.31-0.642x_1+0.032x_2+2.799x_3\\ +1.698x_1^2+5.377x_2^2+4.895x_3^2+5.543x_1{x}_2\\ +1.893x_1{x}_3-0.686x_2{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{2\left(6h\right)}& =5.74-1.195x_1+1.609x_2-5.458x_3\\ +7.112x_1^2+0.037x_2^2+9.608x_3^2+11.9x_1{x}_2\\ -4.042x_1{x}_3+0.98x_2{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{3\left(8h\right)}& =11.548-6.216x_1+3.632x_2-0.354x_3\\ +2.053x_1^2+2.293x_2^2+2.581x_3^2-5.282x_1{x}_2\\ +2.575x_1{x}_3-5.902x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{1\left(1h\right)}& =0.063-0.005x_1-0.001x_2+0.012x_3\\ +0.007x_1^2+0.008x_2^2+0.014x_3^2\\ +0.021x_1{x}_2+0.008x_1{x}_3-0.003x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{2\left(6h\right)}& =0.04-0.002x_1+0.007x_2-0.003x_3\\ +0.006x_1^2-0.004x_2^2+0.009x_3^2\\ +0.013x_1{x}_2-0.008x_1{x}_3+0.004x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{3\left(8h\right)}& =0.039-0.009x_1+0.008x_2+0.002x_3\\ -0.0002x_1^2+0.004x_2^2+0.006x_3^2-0.007x_1{x}_2\\ +0.002x_1{x}_3-0.006x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(s_{12}\right)}_{\left(1h-6h\right)}& =1.89+2.507x_1-0.799x_2+0.299x_3\\ +0.677x_1^2-2.227x_2^2-4.571x_3^2-2.594x_1{x}_2\\ -0.655x_1{x}_3-2.289x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(s_{13}\right)}_{\left(1h-8h\right)}& =1.091+1.603x_{1}1.572x_2+3.023x_3\\ -3.872x_1^2-3.299x_2^2+3.353x_3^2-1.879x_1{x}_2\\ +0.966x_1{x}_3+2.559x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(s_{23}\right)}_{\left(6h-8h\right)}& =-2.945-1.711x_1-1.729x_2-3.296x_3\\ +2.541x_1^2+2.411x_2^2-2.738x_3^2-0.22x_1{x}_2\\ +3.237x_1{x}_3+3.732x_2{x}_3\end{array}$$

Appendix 2

$$\begin{array}{ll}{\mu }_{1\left(0.5h\right)}& =4.844-0.039x_1+0.023x_3-0.006x_7-0.005x_8\\ -0.001x_{10}+0.0001x_1^2-0.00007x_3^2+0.00006x_7^2\\ +0.00002x_8^2+0.00003x_{10}^2+0.0006x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{1\left(0.5h\right)}& =0.71-0.008x_1+0.0001x_3-0.00078x_7\\ +0.006x_8-0.006x_{10}+0.00003x_1^2\\ +0.000006x_3^2+0.00003x_7^2-0.00002x_8^2\\ +0.00004x_{10}^2-0.00003x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{2\left(1h\right)}& =7.644-0.027x_1+0.015x_3-0.01x_7+0.017x_8\\ -0.014x_{10}+0.0001x_1^2+0.000001x_3^2+0.0001x_7^2\\ +0.0002x_8^2+0.0001x_{10}^2+0.0004x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{2\left(1h\right)}& =1.103-0.041x_1-0.027x_3-0.002x_7\\ +0.021x_8+0.006x_{10}+0.0001x_1^2\\ +0.00008x_3^2+0.00002x_7^2-0.00007x_8^2\\ -0.00002x_{10}^2+0.0009x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{3\left(1.5h\right)}& =7.228+0.109x_1+0.018x_3-0.029x_7+0.033x_8\\ -0.035x_{10}-0.0003x_1^2-0.0005x_3^2+0.0003x_7^2\\ +0.0003x_8^2+0.0002x_{10}^2-0.0044x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{3\left(1.5h\right)}& =0.292+0.021x_1+0.035x_3-0.031x_7\\ +0.033x_8-0.004x_{10}-0.00005x_1^2-0.000009x_3^2\\ +0.0002x_7^2-0.0001x_8^2+0.00003x_{10}^2-0.0009x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{4\left(2h\right)}& =8.611+0.165x_1+0.248x_3-0.074x_7\\ +0.074x_8-0.05x_{10}-0.0005x_1^2-0.0007x_3^2\\ +0.0006x_7^2+0.0002x_8^2+0.0003x_{10}^2-0.006x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{4\left(2h\right)}& =1.582-0.082x_1-0.05x_3-0.033x_7+0.058x_8\\ +0.027x_{10}+0.0003x_1^2+0.0002x_3^2\\ +0.0002x_7^2-0.0002x_8^2-0.0001x_{10}^2+0.002x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{5\left(3h\right)}& =12.428+0.207x_1+0.309x_3-0.09x_7\\ +0.089x_8-0.049x_{10}-0.0006x_1^2-0.0008x_3^2\\ +0.0007x_7^2+0.0003x_8^2+0.0004x_{10}^2-0.008x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{5\left(3h\right)}& =1.69-0.078x_1-0.033x_3-0.021x_7+0.052x_8\\ +0.033x_{10}+0.0003x_1^2+0.0001x_3^2\\ +0.0001x_7^2-0.0002x_8^2-0.0001x_{10}^2+0.001x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{6\left(4h\right)}& =16.417+0.287x_1+0.388x_3-0.11x_7\\ +0.126x_8-0.07x_{10}-0.0008x_1^2-0.001x_3^2\\ +0.0009x_7^2+0.0003x_8^2+0.0005x_{10}^2-0.011x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{6\left(4h\right)}& =3.123-0.134x_1-0.074x_3-0.035x_7+0.061x_8\\ +0.053x_{10}+0.0005x_1^2+0.0002x_3^2\\ +0.0002x_7^2-0.0002x_8^2-0.0002x_{10}^2+0.003x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{7\left(6h\right)}& =21.874+0.563x_1+0.691x_3-0.174x_7\\ +0.109x_8-0.084x_{10}-0.002x_1^2-0.002x_3^2\\ +0.001x_7^2+0.0006x_8^2+0.0007x_{10}^2-0.02x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{7\left(6h\right)}& =4.719-0.22x_1-0.104x_3-0.056x_7+0.073x_8\\ +0.105x_{10}+0.0008x_1^2+0.0003x_3^2\\ +0.0003x_7^2-0.0002x_8^2-0.0004x_{10}^2+0.005x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{8\left(8h\right)}& =28.588+0.811x_1+0.963x_3-0.221x_7\\ +0.073x_8-0.11x_{10}-0.002x_1^2-0.003x_3^2\\ +0.001x_7^2+0.0007x_8^2+0.001x_{10}^2-0.03x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{8\left(8h\right)}& =5.417-0.226x_1-0.064x_3-0.072x_7+0.061x_8\\ +0.158x_{10}+0.0008x_1^2+0.0001x_3^2\\ +0.0004x_7^2-0.0002x_8^2-0.0006x_{10}^2+0.004x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{9\left(10h\right)}& =37.1+0.886x_1+1.086x_3-0.249x_7\\ +0.058x_8-0.094x_{10}-0.003x_1^2-0.003x_3^2\\ +0.002x_7^2+0.001x_8^2+0.001x_{10}^2-0.032x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{9\left(10h\right)}& =7.351-0.28x_1-0.085x_3-0.088x_7+0.046x_8\\ +0.201x_{10}+0.001x_1^2+0.0002x_3^2\\ +0.0005x_7^2-0.0002x_8^2-0.0008x_{10}^2+0.005x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{10\left(12h\right)}& =44.362+1.017x_1+1.237x_3-0.229x_7\\ +0.055x_8-0.144x_{10}-0.003x_1^2-0.004x_3^2\\ +0.001x_7^2+0.0006x_8^2+0.001x_{10}^2-0.036x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{10\left(12h\right)}& =7.482-0.267x_1-0.049x_3-0.095x_7+0.055x_8\\ +0.217x_{10}+0.001x_1^2+0.00001x_3^2\\ +0.0005x_7^2-0.0002x_8^2-0.001x_{10}^2+0.004x_1{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{11\left(24h\right)}& =82.688+0.577x_1+0.705x_3-0.056x_7\\ +0.06x_8+0.044x_{10}-0.002x_1^2-0.002x_3^2\\ +0.004x_7^2-0.00004x_8^2+0.0001x_{10}^2-0.02x_1{x}_3\end{array}$$

$$\begin{array}{ll}{V}_{11\left(24h\right)}& =7.503-0.104x_1-0.025x_3-0.097x_7-0.005x_8\\ -0.004x_{10}+0.0005x_1^2+0.00004x_3^2\\ +0.0006x_7^2-0.0001x_8^2-0.0001x_{10}^2+0.001x_1{x}_3\end{array}$$

Appendix 3

$$\begin{array}{ll}{\mu }_{1\left(2h\right)}& =12.986-2.16x_1-x_2+0.68x_3+0.121x_1^2-0.279x_2^2\\ +0.221x_3^2+0.038x_1{x}_2+0.038x_1{x}_3+0.163x_2{x}_3\end{array}$$

$$\begin{array}{ll}{v}_{1\left(2h\right)}& =1.274+0.057x_1+0.33x_2-0.235x_3-0.064x_1^2\\ +0.056x_2^2-0.298x_3^2-0.002x_1{x}_2+0.426x_1{x}_3\\ +0.292x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{1\left(2h\right)}& =0.082+0.013x_1+0.017x_2-0.015x_3\\ -0.006x_1^2+0.011x_2^2-0.016x_3^2+0.002x_1{x}_2\\ +0.014x_1{x}_3+0.015x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{2\left(4h\right)}& =25.121-5.2x_1-2x_2+1.43x_3+0.47x_1^2-0.331x_2^2\\ +0.619x_3^2+0.163x_1{x}_2+0.063x_1{x}_3+0.338x_2{x}_3\end{array}$$

$$\begin{array}{ll}{v}_{2\left(4h\right)}& =1.747+0.112x_1+0.004x_2-0.564x_3-1.017x_1^2\\ +1.813x_2^2-0.732x_3^2-0.442x_1{x}_2\\ +0.185x_1{x}_3-0.185x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{2\left(4h\right)}& =0.046+0.01x_1+0.002x_2-0.011x_3\\ -0.013x_1^2+0.025x_2^2-0.009x_3^2-0.004x_1{x}_2\\ +0.0004x_1{x}_3-0.004x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{3\left(6h\right)}& =42.938-7.27x_1-2.87x_2+2.31x_3\\ -0.257x_1^2-0.257x_2^2-0.057x_3^2+0.913x_1{x}_2+0.463x_1{x}_3\\ -0.688x_2{x}_3\end{array}$$

$$\begin{array}{ll}{v}_{3\left(6h\right)}& =3.412+0.072x_1+0.965x_2+0.052x_3\\ +3.869x_1^2-1.106x_2^2-1.351x_3^2-1.126x_1{x}_2\\ +0.936x_1{x}_3-0.049x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{3\left(6h\right)}& =0.042+0.009x_1+0.008x_2-0.002x_3\\ -0.023x_1^2-0.005x_2^2-0.007x_3^2-0.006x_1{x}_2\\ +0.004x_1{x}_3+0.001x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{4\left(9h\right)}& =67.278-11.37x_1-3.02x_2+3.27x_3-2.541x_1^2\\ +3.709x_2^2-3.841x_3^2+0.125x_1{x}_2+0.825x_1{x}_3\\ +0.05x_2{x}_3\end{array}$$

$$\begin{array}{ll}{v}_{4\left(9h\right)}& =3.563+0.311x_1-0.064x_2+0.085x_3-0.895x_1^2\\ +0.32x_2^2+0.425x_3^2+0.523x_1{x}_2-0.208x_1{x}_3\\ -0.09x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{4\left(9h\right)}& =0.027+0.007x_1+0.002x_2-0.002x_3\\ -0.001x_1^2+0.001x_2^2+0.003x_3^2+0.003x_1{x}_2\\ -0.002x_1{x}_3-0.001x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{5\left(12h\right)}& =82.395-12.84x_1-5.25x_2+3.8x_3-0.567x_1^2\\ -0.417x_2^2+0.333x_3^2-0.675x_1{x}_2+0.625x_1{x}_3\\ +0.125x_2{x}_3\end{array}$$

$$\begin{array}{ll}{v}_{5\left(12h\right)}& =3.944-0.428x_1+0.038x_2-0.142x_3+1.018x_1^2\\ -1.592x_2^2-0.662x_3^2+0.705x_1{x}_2-0.065x_1{x}_3\\ +0.643x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{5\left(12h\right)}& =0.024+0.002x_1+0.002x_2-0.002x_3\\ +0.004x_1^2-0.006x_2^2-0.002x_3^2+0.003x_1{x}_2\\ -0.001x_1{x}_3+0.003x_2{x}_3\end{array}$$

Appendix 4

$$\begin{array}{ll}{\mu }_{1\left(0.5h\right)}& =31.153-3.546x_1-3.884x_2+3.243x_3\\ +0.667x_1^2+1.874x_2^2-3.391x_3^2\\ +2.897x_1{x}_2-0.767x_1{x}_3+1.175x_2{x}_3\end{array}$$

$$\begin{array}{ll}{v}_{1\left(0.5h\right)}& =0.669-0.456x_1-0.45x_2-0.839x_3+1.542x_1^2\\ -1.429x_2^2+2.026x_3^2-1.309x_1{x}_2-1.167x_1{x}_3\\ +0.649x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{1\left(0.5h\right)}& =0.028+0.002x_1+0.0004x_2-0.016x_3\\ +0.01x_1^2-0.01x_2^2+0.022x_3^2-0.016x_1{x}_2\\ -0.01x_1{x}_3+0.002x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{2\left(3.5h\right)}& =64.474-6.603x_1-4.648x_2+3.1x_3-0.977x_1^2\\ +4.658x_2^2+1.287x_3^2-1.168x_1{x}_2-0.65x_1{x}_3\\ -0.705x_2{x}_3\end{array}$$

$$\begin{array}{ll}{v}_{2\left(3.5h\right)}& =0.841-0.063x_1+0.215x_2+0.12x_3-0.173x_1^2\\ +0.765x_2^2+0.048x_3^2-0.084x_1{x}_2-0.56x_1{x}_3\\ -0.371x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{2\left(3.5h\right)}& =0.011+0.001x_1+0.003x_2-0.001x_3\\ -0.001x_1^2+0.007x_2^2-0.0003x_3^2-0.00003x_1{x}_2\\ -0.004x_1{x}_3-0.002x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\mu }_{3\left(8h\right)}& =92.466-4.383x_1-2.878x_2+1.811x_3-1.242x_1^2\\ +2.206x_2^2-0.987x_3^2-1.1x_1{x}_2+0.168x_1{x}_3\\ +2.018x_2{x}_3\end{array}$$

$$\begin{array}{ll}{v}_{3\left(8h\right)}& =0.895-0.192x_1+0.213x_2-0.302x_3\\ +0.029x_1^2-0.564x_2^2+0.786x_3^2-0.135x_1{x}_2\\ +0.088x_1{x}_3-0.284x_2{x}_3\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{3\left(8h\right)}& =0.01-0.001x_1+0.001x_2-0.001x_3\\ +0.0001x_1^2-0.004x_2^2+0.004x_3^2-0.001x_1{x}_2\\ +0.001x_1{x}_3-0.001x_2{x}_3\end{array}$$

Appendix 5

$$\begin{array}{ll}{\mu }_{1\left(1h\right)}& =37.191-7.918x_1-3.955x_2+1.148x_1^2-1.432x_2^2\\ -0.558x_1{x}_2\end{array}$$

$$\begin{array}{ll}{v}_{1\left(1h\right)}& =1.957+0.862x_1-0.693x_2-0.105x_1^2-0.04x_2^2\\ -1.32x_1{x}_2\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{1\left(1h\right)}& =0.038+0.015x_1-0.003x_2-0.002x_1^2-0.001x_2^2\\ -0.012x_1{x}_2\end{array}$$

$$\begin{array}{ll}{\mu }_{2\left(5h\right)}& =75.29-6.358x_1-8.795x_2+1.035x_1^2-1.345x_2^2\\ +0.745x_1{x}_2\end{array}$$

$$\begin{array}{ll}{v}_{2\left(5h\right)}& =5.129+0.25x_1+0.915x_2-2.223x_1^2-0.583x_2^2\\ -1.18x_1{x}_2\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{2\left(5h\right)}& =0.031+0.003x_1+0.006x_2-0.009x_1^2\\ -0.002x_2^2-0.005x_1{x}_2\end{array}$$

$$\begin{array}{ll}{\mu }_{3\left(10h\right)}& =89.216-8.49x_1-7.528x_2+3.797x_1^2-1.728x_2^2\\ -3.195x_1{x}_2\end{array}$$

$$\begin{array}{ll}{v}_{3\left(10h\right)}& =3.026-0.145x_1-1.292x_2+2.372x_1^2-1.968x_2^2\\ +0.75x_1{x}_2\end{array}$$

$$\begin{array}{ll}{\left(\frac{s}{m}\right)}_{3\left(10h\right)}& =0.017+0.002x_1-0.003x_2+0.007x_1^2\\ -0.004x_2^2-0.003x_1{x}_2\end{array}$$

$$\begin{array}{ll}{f}_2\mathit{value}& =50.157+7.52x_1+9.473x_2-5.26x_1^2-1.49x_2^2\\ -0.66x_1{x}_2\end{array}$$

Appendix 6

$$\begin{array}{ll}{\mu }_{1\left(1h\right)}& =20.778-3.317x_1-4.017x_2+0.183x_1^2-0.917x_2^2\\ -0.325x_1{x}_2\end{array}$$

$$\begin{array}{ll}{\mu }_{2\left(4h\right)}& =38.678-4.5x_1-5.7x_2+1.583x_1^2-1.467x_2^2\\ -1.425x_1{x}_2\end{array}$$

$$\begin{array}{ll}{\mu }_{3\left(12h\right)}& =68.822-5.483x_1-5.5x_2+2.317x_1^2-1.333x_2^2\\ +0.15x_1{x}_2\end{array}$$

$${\mu }_{4\left(t_{50}\right)}=6.222+x_1+1.167x_2-0.333x_1^2+0.167x_2^2$$

$$\begin{array}{ll}{\mu }_{5\left(\mathit{MDT}\right)}& =8.222+0.767x_1+0.933x_2-0.333x_1^2\\ +0.267x_2^2-0.1x_1{x}_2\end{array}$$

$$\begin{array}{ll}{\mu }_{6\left(f2\right)}& =68.556+11.183x_1+11.45x_2-2.483x_1^2\\ -3.583x_2^2-1.525x_1{x}_2\end{array}$$