Central Composite Design, quadratic model fitting, and desirability-based optimization for Analytical Quality by Design
Every response surface study begins with a design decision that matters more than any statistical technique that follows: which factors to vary, over what ranges, and which response actually reflects process performance. Get this wrong and no amount of sophisticated regression will rescue the study. Design of Experiments (DoE) exists precisely because one-factor-at-a-time (OFAT) tinkering — the default instinct of most experimentalists — systematically fails to characterize how real processes behave.
OFAT experimentation holds every factor fixed except one, sweeps that factor across its range, records the optimum, fixes it there, and moves to the next factor. It feels intuitive and is taught informally in most laboratories — and it is quietly wrong for two structural reasons.
First, OFAT cannot detect interactions. If the true optimum temperature depends on which pH you are running at (a very common situation in reaction chemistry and bioprocessing), then optimizing temperature at one pH and pH at some other temperature will converge on a point that is optimal for neither combination. The experimenter has no way to even notice the problem — the data was never collected in a way that could reveal it.
Second, OFAT wastes information. Each run in a one-factor-at-a-time sweep informs you about one factor only. Each run in a factorial design informs you simultaneously about every factor and every interaction, because all factors change together in a structured, balanced pattern. A well-built response surface design routinely extracts more usable information from 13–20 runs than a comparable OFAT campaign extracts from 40–60.
Studies comparing the two approaches on real processes consistently find that OFAT-derived "optima" leave 10–30% of achievable response on the table relative to the true optimum recovered by response surface methodology (RSM) — because the true peak of the response surface almost never sits on the ridge line that OFAT happens to trace out.
A response surface study starts with risk assessment: which process inputs plausibly affect the response, and which of those are continuous and controllable enough to manipulate deliberately? For a typical reaction or bioprocess optimization, three candidates are common — Temperature, pH, and Reaction Time — but visualizing a surface in more than two dimensions is difficult, so many studies fix the least influential factor (often Time) at its center level and build a 2-factor CCD in Temperature and pH, exactly as this simulation does.
Each factor is assigned a low level (−1), a center level (0), and a high level (+1) in coded units, plus extended axial levels (±α) that extrapolate slightly beyond the cube to estimate curvature. Coding factors this way — (actual − center) / (half-range) — makes the regression numerically well-conditioned and lets coefficients from factors measured in wildly different units (°C vs. pH units) be compared directly on the same scale.
The response must be a single, precisely measurable quantity that reflects what "good" means for the process — Yield (%), Purity, potency, or a critical quality attribute in a QbD context. When more than one response matters (e.g., Yield and Impurity simultaneously) they are each modeled separately and reconciled later with a desirability function (Stage 5).
Coded units are a convention, not a cosmetic detail: −1.414 and +1.414 (α = 2^(k/4) for a rotatable CCD with k=2 factors) are the axial points that extend just beyond the ±1 cube along each axis alone, letting the design estimate the pure quadratic curvature term for each factor independently of the corner points.
A response surface cannot be characterized by a straight line — it needs at least three levels per factor to detect curvature. The Central Composite Design (CCD), introduced by Box and Wilson in 1951, is the workhorse geometry for this task: it augments a two-level factorial cube with axial "star" points and center-point replicates to estimate a full quadratic model with remarkably few runs.
A CCD is built from three interleaved point sets, each with a distinct statistical job:
• Factorial (corner) points — the 2^k vertices of the design cube (4 points for k=2: (−1,−1), (1,−1), (−1,1), (1,1)). These estimate the linear main effects and the two-factor interaction term with maximum efficiency, exactly as in a classical two-level factorial.
• Axial / star points — 2k points sitting on each factor axis at ±α with all other factors at 0 (for k=2: (−1.414,0), (1.414,0), (0,−1.414), (0,1.414)). These are what make curvature estimable: without them, the design could only fit a planar (linear + interaction) model, never a bowl or a saddle.
• Center points — several replicate runs at the exact center (0,0) of the design space. These serve two purposes: they provide a model-independent estimate of pure experimental error (needed for the lack-of-fit test in Stage 4), and their number can be tuned to balance the design's prediction variance across the whole region — more replicates sharpen the error estimate but add runs.
For two factors this yields the classic 13-run CCD (4 + 4 + 5) that appears in almost every RSM textbook, and which the "Center Point Replicates" slider in this simulation lets you resize between 3 and 8 replicates.
When the ±1.414 axial extremes of a CCD are impractical or unsafe (a temperature or pressure combination that risks damaging equipment or degrading product), the Box–Behnken design (BBD) offers an alternative geometry. Instead of a cube-plus-star, BBD places points at the midpoints of the edges of the design cube, combined with center points — every factor only ever takes 3 distinct levels (−1, 0, +1), and no run ever sets every factor to its extreme simultaneously.
BBD is typically more run-efficient than a CCD for three or more factors and is the preferred choice in pharmaceutical and bioprocess development, where the corners of a CCD cube can represent combinations (e.g., maximum temperature at minimum pH) that were never physically tested and may be outside a safe or even physically meaningful operating envelope.
| Product | Indication | Trial Design | Key Result |
|---|---|---|---|
| OFAT sweep | 1 factor at a time | Sequential univariate sweeps; each run informs only one factor | Simple, but blind to interactions and curvature |
| Full 3-level factorial | All factors, 3 levels | 3^k exhaustive grid; captures curvature completely | Complete but expensive — 27 runs for k=3 |
| Central Composite (CCD) | Cube + axial + center | 2^k + 2k + nc runs; rotatable, efficient quadratic fit | Gold standard for RSM; minimal runs for full quadratic |
| Box–Behnken (BBD) | Cube edge midpoints + center | No extreme corner combinations; 3 levels per factor | Safer operating envelope, run-efficient for k≥3 |
A design matrix is only a plan until it is executed. Turning 13 rows of coded factor levels into 13 measured Yield values requires disciplined execution: randomized run order, awareness of nuisance variables, and honest measurement — because every downstream statistic (R², lack-of-fit, the fitted surface itself) inherits whatever noise or bias creeps in at this stage.
If the 13 runs of a CCD were executed in their tabulated order — all factorial points first, then all axial points, then all center points — any drift over time (a reagent degrading, a sensor calibration slipping, an operator learning the procedure) would be perfectly confounded with the point type, corrupting exactly the curvature estimate the design was built to capture.
Randomizing the execution order breaks this link: nuisance drift is scattered randomly across corner, [[-ALPHA,0],[ALPHA,0],[0,-ALPHA],[0,ALPHA]], and center points instead of being concentrated on one type, so it inflates residual noise slightly rather than systematically biasing specific model terms. When a known nuisance factor cannot be randomized away (e.g., only one batch of raw material fits in a single day), it is instead handled by blocking — deliberately grouping runs so the nuisance factor becomes a controlled block effect in the model rather than an unmodeled confound.
Center points are typically interspersed throughout the run order rather than run consecutively. Repeating them at intervals throughout the campaign lets the experimenter detect a time-trend directly — if early and late center-point replicates disagree noticeably, something in the process drifted during the study.
As each run completes, its measured response is entered into the design matrix alongside its coded factor levels — the tabular structure that both the corner-point and axial-point rows share is what makes the least-squares fit in Stage 4 a routine matrix computation rather than a bespoke calculation per experiment.
Before fitting any model, the completed matrix is screened for obvious problems: transcription errors, assay outliers, and runs that silently violated their intended factor levels (a furnace that overshot its temperature setpoint, for instance). Center-point replicates are particularly useful here — because they are nominally identical runs, their spread gives an immediate, model-free sense of how much random noise to expect before the regression is even fit.
With a completed design matrix in hand, response surface methodology fits a second-order (quadratic) polynomial relating the response to the coded factors. Unlike a simple linear regression, this model captures three distinct kinds of process behavior at once: how much each factor matters on its own, how factors modify each other's effect, and whether the response curves rather than climbing in a straight line.
For two coded factors X1 and X2, the full quadratic response surface model is:
Y = b0 + b1·X1 + b2·X2 + b12·X1X2 + b11·X1² + b22·X2² + ε
• b0 — the intercept, the predicted response at the center point • b1, b2 — linear main-effect coefficients: how steeply Y responds to each factor near the center • b12 — the two-factor interaction coefficient: how much the effect of X1 depends on the current level of X2 (and vice versa) — the term OFAT experimentation can never estimate • b11, b22 — pure quadratic (curvature) coefficients: whether the response bends into a peak, a valley, or a saddle as each factor moves away from center — the terms that require axial points to estimate • ε — residual error, everything the model does not explain
Coefficients are estimated by ordinary least squares: b = (XᵀX)⁻¹XᵀY, where X is the design matrix expanded with interaction and squared columns. Because a CCD is built to be (nearly) orthogonal, the coefficient estimates are largely uncorrelated with one another — each term can be interpreted with minimal risk of one estimate distorting another.
A fitted model is not trustworthy until it has passed three checks:
• R² and adjusted R² — the fraction of response variance explained by the model. R² alone can be inflated by adding terms indiscriminately; adjusted R² penalizes unnecessary complexity and is the more honest number to report.
• Lack-of-fit F-test — this is where replicated center points earn their keep. Total residual error is decomposed into "pure error" (the scatter among identical center-point replicates — measurement noise the model could never explain) and "lack-of-fit" (systematic deviation between the model and the data beyond that noise floor). A non-significant lack-of-fit p-value (> 0.05) means the quadratic form is adequate; a significant one means the true surface has structure — cubic terms, a discontinuity — that a quadratic cannot capture.
• Residual analysis — plotting residuals against predicted values, against run order, and on a normal probability plot. A random, structureless scatter supports the model's assumptions (constant variance, independence, normality); a funnel, curve, or trend in any of these plots signals a violated assumption that undermines the model's predictions and confidence intervals.
A model can have an excellent R² and still fail the lack-of-fit test — high R² only says the model explains most of the variance actually present in the data, while lack-of-fit specifically asks whether that unexplained residual is pure noise or a symptom of missing curvature. Both checks are needed before trusting the surface for optimization.
A validated quadratic model is a compact, continuous description of the entire factor space — not just the 13 points that were actually run. Contour and 3-D surface plots make that description visually interpretable, and a stationary-point analysis or desirability function turns it into an actionable recommendation: the specific factor combination to run next.
A contour plot slices the fitted surface into bands of equal predicted response — concentric rings around a peak indicate a true maximum; elongated, parallel bands indicate a ridge where many factor combinations give nearly the same response (useful for finding a robust operating region, not just a single fragile optimum); a saddle-shaped pattern (rings curving oppositely along the two axes) signals that the "stationary point" found by calculus is neither a maximum nor a minimum but a pass between two opposing slopes.
Formally, the nature of the stationary point is determined by the eigenvalues of the matrix of quadratic coefficients (B, built from b11, b22, and b12/2): both eigenvalues negative means a true maximum, both positive a true minimum, and mixed signs a saddle point. This eigenvalue analysis (canonical analysis) is what separates a rigorous RSM optimum from simply eyeballing the highest contour band.
Real processes are rarely optimized on a single response. Maximizing Yield while simultaneously minimizing an Impurity and holding Purity above a specification requires reconciling several fitted surfaces at once — exactly the problem the Derringer–Suich desirability function solves.
Each response Yi is transformed into an individual desirability di ∈ [0,1]: di = 1 where Yi is ideal, di = 0 where Yi is unacceptable, with a smooth ramp (often shaped by a weighting exponent) in between. The individual desirabilities are then combined into a single overall desirability D = (d1 × d2 × ... × dn)^(1/n) — a geometric mean, chosen specifically because it collapses to zero if any single response is unacceptable, no matter how good the others are. The factor combination that maximizes D is the recommended joint optimum across every response simultaneously.
Response surface methodology is not just an optimization trick — it is the quantitative engine behind Quality by Design (QbD), the ICH Q8(R2) framework now standard in pharmaceutical process development. Once a validated quadratic model links critical process parameters to critical quality attributes, the same surface used to find a single optimum can instead be used to map an entire design space: the multidimensional combination of factor ranges within which the model predicts quality attributes will remain within specification with high confidence, typically visualized as an overlay of multiple contour plots — one per quality attribute — intersected into a single acceptable region.
Operating anywhere inside a regulator-accepted design space is not considered a process change requiring re-approval, which is the core practical payoff of RSM-driven QbD: it converts single-point process optimization into a robust, flexible operating region, verified once and defensible for the life of the product.
ICH Q8(R2) defines design space as "the multidimensional combination and interaction of input variables and process parameters that have been demonstrated to provide assurance of quality." A response surface — corner points, axial points, center-point replicates, quadratic fit, and all — is the experimental and statistical machinery that makes such a claim defensible.