Mixture Design

This tutorial illustrates the use of a mixture design to investigate the formulation of a pharmaceutical tablet.

The tablet formulation contains several fixed ingredients:

• drug substance: \(4.0\%\);

• disintegrant: \(0.2\%\);

• lubricant A: \(4.0\%\);

• lubricant B: \(1.0\%\).

The remaining \(90.8\%\) of the formulation is composed of three possible binding agents:

• Monohydrated :math:`alpha`-lactose, :math:`X_1`;

• Anhydrous :math:`beta`-lactose, :math:`X_2`;

• Modified rice starch, :math:`X_3`.

Only the relative proportions of these three components are investigated. Their normalized fractions must satisfy the mixture constraint:

\[X_1 + X_2 + X_3 = 1\]

Unlike ordinary process factors, mixture components cannot vary independently. Increasing the proportion of one component necessarily reduces the proportion of one or both of the remaining components.

Experimental design

For a three-component mixture, the experimental region is represented by an equilateral triangle:

• the vertices correspond to the pure components;

• the edges correspond to binary mixtures;

• the interior points correspond to ternary mixtures.

The design contains seven experiments:

• three pure-component formulations;

• three binary mixtures containing equal proportions of two components;

• one ternary mixture containing equal proportions of all three components.

Response variables

Two responses are measured for each formulation:

• :math:`Y_1`: breaking load, representing the maximum force the tablet can withstand before breaking;

• :math:`Y_2`: dissolution rate.

The aim is to identify the best compromise between the two responses. A formulation with a high breaking load is desirable because it produces mechanically resistant tablets, while a high dissolution rate is required for rapid drug release.

These objectives may conflict, since formulations with greater mechanical strength may dissolve more slowly.

Model specification

A special cubic mixture model is considered for each response:

\[Y = \beta_1X_1 + \beta_2X_2 + \beta_3X_3 + \beta_{12}X_1X_2 + \beta_{13}X_1X_3 + \beta_{23}X_2X_3 + \beta_{123}X_1X_2X_3 + \varepsilon\]

where:

• \(\beta_1\), \(\beta_2\), and \(\beta_3\) describe the responses associated with the pure components;

• \(\beta_{12}\), \(\beta_{13}\), and \(\beta_{23}\) describe non-additive blending effects in the binary mixtures;

• \(\beta_{123}\) describes the ternary blending effect;

• \(\varepsilon\) represents the experimental error.

No intercept is included because a mixture composition with all components equal to zero is not physically possible under the constant-sum constraint.

The model contains seven coefficients, equal to the number of experimental formulations.

In this notebook, doetools is used to:

1. define the three mixture components;

2. generate the seven-point mixture design;

3. inspect and visualize the triangular experimental region;

4. specify the special cubic mixture model;

5. fit separate models for breaking load and dissolution rate;

6. evaluate the compromise between the two formulation responses.

This example is adapted from Leardi (2009), Experimental design in chemistry: A tutorial, published in Analytica Chimica Acta, DOI: https://doi.org/10.1016/j.aca.2009.06.015.

# This is needed to display plotly figures in the Jupyter notebook, not for the library itself
import plotly.io as pio

pio.renderers.default = "notebook"

# Fundamental imports
from doetools import MixtureFactor, ModelTerms, SimplexLatticeDesign

---

1. Design Initialization and Exploration

1.1 Factor Definition and Design Creation

Define the three mixture components and create the Simplex Lattice Design object.

Each component is a MixtureFactor with lower_bound=0 and upper_bound=1, meaning each can range from 0% to 100% of the blend. The key constraint is that all proportions must always sum to 1:

\[x_1 + x_2 + x_3 = 1, \quad x_i \geq 0\]

The {m, 3} simplex lattice places points at every combination where each component takes values in \(\{0, \frac{1}{m}, \frac{2}{m}, \ldots, 1\}\). With m = 2 we get the vertices (pure blends), edge midpoints (binary blends), and optionally center points (ternary blend), which is enough to fit a second-order Scheffé polynomial.

# Define 3 mixture components (ligands)
factors = {
    "X1" : MixtureFactor(lower_bound=0, upper_bound=1),  # Ligand 1
    "X2" : MixtureFactor(lower_bound=0, upper_bound=1),  # Ligand 2
    "X3" : MixtureFactor(lower_bound=0, upper_bound=1)   # Ligand 3
}

# Create the Simplex Lattice Design with m=2 and 1 center point
design = SimplexLatticeDesign(factors=factors, m=2, center_points=1)

1.2 Design Matrix

Display the complete design matrix showing all experimental runs in real units (component proportions).

design_matrix = design.get_design_matrix()
design_matrix

	X1	X2	X3
0	0.000000	0.000000	1.000000
1	0.000000	0.500000	0.500000
2	0.000000	1.000000	0.000000
3	0.500000	0.000000	0.500000
4	0.500000	0.500000	0.000000
5	1.000000	0.000000	0.000000
6	0.333333	0.333333	0.333333

Visualize the experimental design on a ternary plot (simplex).

fig00 = design.plot_design(ax1="X1", ax2="X2", ax3="X3")
fig00

1.3 Export Design for Experimentation

Export the design matrix to an Excel file for use in the laboratory. The exported file:

• Contains all experimental runs with the required component proportions

• Includes empty columns (Y1, Y2) for recording measured breaking load and dissolution rate

Next step: Prepare the tablet batches following the design matrix, measure breaking load and dissolution rate, and record the values in the Excel file.

design.export_experiments(responses=["Y1", "Y2"],
                          randomize=False,
                          destination="03_mix_design.xlsx")

---

2. Model Specification

Define the regression model structure for mixture response surface analysis. For mixture experiments the conventional choice is the Scheffé polynomial, which respects the mixture constraint \(\sum x_i = 1\) and avoids redundant terms.

The second-order (quadratic) Scheffé model for three components is:

\[Y = \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + \beta_{12} x_1 x_2 + \beta_{13} x_1 x_3 + \beta_{23} x_2 x_3 + \beta_{123} x_1 x_2 x_3\]

Where:

• βᵢ (Linear blending terms): Effect of each pure component — the predicted response when only component i is present

• βᵢⱼ (Binary synergy/antagonism terms): Departure from linear blending — positive values indicate synergistic mixing, negative values indicate antagonism

• Ternary interaction: Complex ternary interactions between components

Note that this model has no intercept term — the intercept is already absorbed into the linear blending coefficients due to the mixture constraint. This is automatically handled by setting intercept=False.

model_terms = ModelTerms(intercept=False,
                         scheffe_pol_order=3)

design.set_model_terms(model_terms)

---

3. Leverage Analysis

3.1 Leverage Values

Leverage (\(h_{ii}\)) quantifies how much each experimental point influences the fitted model. It represents the diagonal elements of the hat matrix: \(\mathbf{H} = \mathbf{X}(\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\).

In this case the design is saturated, so the leverage value of experimental points is one.

leverages = design.get_leverages()
leverages

	Leverage
0	1.0
1	1.0
2	1.0
3	1.0
4	1.0
5	1.0
6	1.0

---

3. Response Analysis and Model Fitting

3.1 Import Experimental Results

After completing the tablet preparation and testing, import the measured response data. The complete analysis workflow:

1. Import responses: Read the Excel file containing breaking load (Y1) and dissolution rate (Y2)

2. Fit MLR model: Estimate the Scheffé polynomial coefficients via least squares

3. Set optimization criteria: Define constraints and objectives for each response

   • Y1 (Breaking load): Maximize, lower limit = 50 — tablets must have sufficient mechanical resistance to survive handling and packaging

   • Y2 (Dissolution rate): Maximize, lower limit = 50 — adequate drug release must be ensured within 30 minutes

# Import Responses and Fit MLR Model
design.import_responses(source="03_mix_responses.xlsx")
design.compute_mlr_model()

# Set Response Conditions
design.set_response_conditions(lower_limits=[50, 50],
                               upper_limits=[None, None],
                               maximize=[True, True])

# Export the Mixture Deisign Report to PDF
design.get_model_summary_pdf("03_mix_report.pdf")

5.2 Regression Coefficients — Y1 (Breaking Load)

design.get_coefficient_summary("Y1")

	Variable	Coefficient	Std_Error	Conf_Int	Lower_CI	Upper_CI	p_value
0	X1	31.0	NaN	NaN	NaN	NaN	NaN
1	X2	113.0	NaN	NaN	NaN	NaN	NaN
2	X3	38.0	NaN	NaN	NaN	NaN	NaN
3	X1:X2	-120.0	NaN	NaN	NaN	NaN	NaN
4	X1:X3	18.0	NaN	NaN	NaN	NaN	NaN
5	X2:X3	-22.0	NaN	NaN	NaN	NaN	NaN
6	X1:X2:X3	354.0	NaN	NaN	NaN	NaN	NaN

5.3 Regression Coefficients — Y2 (Dissolution Rate)

design.get_coefficient_summary("Y2")

	Variable	Coefficient	Std_Error	Conf_Int	Lower_CI	Upper_CI	p_value
0	X1	44.0	NaN	NaN	NaN	NaN	NaN
1	X2	26.0	NaN	NaN	NaN	NaN	NaN
2	X3	52.0	NaN	NaN	NaN	NaN	NaN
3	X1:X2	24.0	NaN	NaN	NaN	NaN	NaN
4	X1:X3	76.0	NaN	NaN	NaN	NaN	NaN
5	X2:X3	-36.0	NaN	NaN	NaN	NaN	NaN
6	X1:X2:X3	195.0	NaN	NaN	NaN	NaN	NaN

---

6. Response Surface Analysis

6.1 Response Surface — Y1 (Breaking Load)

Visualize the fitted Scheffé model for Y1 as a ternary contour map and 3D surface. The ternary plot is the natural representation for mixture data.

y1_cp, y1_sur = design.plot_response(ax1="X1", ax2="X2", ax3="X3",
                                     response="Y1",
                                     feasible_region=True)

Ternary Contour Plot — Y1

y1_cp

6.2 Response Surface — Y2 (Dissolution Rate)

Response surface for Y2 (dissolution rate). Compare the topology with Y1 to identify trade-offs between breaking load and dissolution rate. Regions where both responses are high define the optimal formulation space.

y2_cp, y2_sur = design.plot_response(ax1="X1", ax2="X2", ax3="X3",
                                     response="Y2",
                                     feasible_region=True)

Ternary Contour Plot — Y2

y2_cp

6.3 Multi-Objective Optimization

The critical plot for formulation decision-making! This overlay simultaneously displays both responses to identify the feasible formulation region where all constraints are met.

How to read this plot:

• Highlighted/colored region: Compositions satisfying BOTH Y1 ≥ 50 AND Y2 ≥ 50

Decision strategy:

• If a feasible region exists → Select a formulation within this zone; prefer compositions near the center of the feasible region for robustness against raw material variability

• If no feasible region exists → Relax constraints, explore extended component bounds, or consider adding a fourth component

For this tablet formulation: choose the X1/X2/X3 proportions within the feasible zone that best balance manufacturing practicality with performance.

y3_cp, y3_sur = design.plot_response(ax1="X1", ax2="X2", ax3="X3",
                                     response="Y1",
                                     second_response="Y2",
                                     feasible_region=True)

Ternary Overlay Plot

y3_cp

6.4. Mixture Trace

# The mixture trace plot shows how the predicted response changes as each component
# is varied from the reference formulation, while the remaining components are adjusted
# to satisfy the mixture constraint.

design.plot_mixture_trace("Y1")