Morris Screening for PBPK Parameters

Why Use Morris Screening?

A PBPK model may contain many uncertain physiological and drug-specific parameters, while the outcome of interest may be a single exposure metric, such as area under the concentration-time curve (AUC). A local derivative can describe how the output responds near one parameter vector, but it may miss changes elsewhere in the parameter ranges or changes in effect direction. Morris screening addresses this by evaluating parameter changes at multiple locations in the input space. It is designed to help distinguish parameters that appear influential over specified ranges from those that appear less influential, at a lower computational cost than evaluating every possible combination.

The method is a screening tool, not a complete account of global output variation. Its results depend on the selected parameter ranges, the output being examined, the trajectory sample, and the numerical reliability of the PBPK simulations. A parameter can be influential for one endpoint and unimportant for another.

Definition: For a scalar model output $f(\mathbf{z})$, where each $z_i$ is a normalized input ranging from 0 to 1, the elementary effect for input $i$ at a point $\mathbf{z}$ is $$ EE_i(\mathbf{z})=\frac{f(\mathbf{z}+\Delta e_i)-f(\mathbf{z})}{\Delta}, $$ provided both points are in the allowed input grid. Here $e_i$ is the unit vector in coordinate $i$, and $\Delta$ is the chosen grid step. The elementary effect is a finite change in output per unit change in normalized input.

Set Up the Parameter Space and Trajectories

Before calculating effects, specify the scalar PBPK output and define a range for every parameter to be screened. For a parameter $\theta_i$ with lower and upper bounds $L_i$ and $U_i$, a linear scaling is $z_i=(\theta_i-L_i)/(U_i-L_i)$. Thus $z_i=0$ corresponds to $L_i$, and $z_i=1$ corresponds to $U_i$. If a parameter is more naturally explored on a logarithmic scale, define and document a scaling appropriate to that scale before constructing trajectories. In either case, the normalized coordinates must represent the intended range of plausible values.

A common Morris design uses $p$ equally spaced levels from 0 to 1, with even $p$, and a step size $\Delta=p/[2(p-1)]$. For example, with $p=4$, the grid levels are 0, $1/3$, $2/3$, and 1, and $\Delta=2/3$. The starting coordinates are restricted so that a step of size $\Delta$ in the selected direction remains on the grid and within the range. A trajectory changes one coordinate at each step: first choose a starting point, then a random ordering of the $k$ parameters, and then a valid direction for each change. A trajectory therefore requires $k+1$ model evaluations for $k$ parameters.

Repeating the construction with different starting points, orderings, and directions gives $r$ elementary effects for each parameter. All parameters are varied one at a time along any single trajectory, but effects are sampled at different locations across trajectories. This is why Morris screening can reveal location-dependent behavior that a single one-at-a-time calculation at a reference vector does not.

Choose the target.
State the PBPK model output, such as AUC in mg·h/L or a dimensionless exposure ratio. Analyze each output separately.

Define and scale ranges.
Specify plausible bounds and map each parameter to a normalized coordinate. Keep the ranges and scaling fixed for the analysis.

Generate trajectories.
Sample valid grid starting points, coordinate orderings, and directions; evaluate the model at every point on each trajectory.

Calculate and summarize effects.
Compute one finite-difference effect per parameter per trajectory, then examine its signed mean, mean absolute value, and variability.

Summaries and Their Interpretation

For parameter $i$, let $EE_{ij}$ be its elementary effect from trajectory $j$, for $j=1,\ldots,r$. Three summaries are commonly used. The signed mean $\mu_i$ can show the average direction of change, but positive and negative effects can cancel. The mean absolute effect $\mu_i^*$ avoids that cancellation and indicates the average magnitude of the response over the sampled locations. The sample standard deviation $\sigma_i$ describes how much the effects vary across trajectories.

\mu_i=\frac{1}{r}\sum_{j=1}^{r}EE_{ij},\qquad \mu_i^*=\frac{1}{r}\sum_{j=1}^{r}|EE_{ij}|,\qquad \sigma_i=\sqrt{\frac{1}{r-1}\sum_{j=1}^{r}(EE_{ij}-\mu_i)^2}.

Because inputs are normalized, an elementary effect is expressed in units of the output per unit change across the normalized input scale. If the output is AUC in mg·h/L, then each effect and its summaries are also in mg·h/L per normalized input unit. A large $\mu_i^*$ suggests that changing parameter $i$ across the specified range often changes the chosen output substantially. A small value suggests relatively little effect over those ranges and sampled locations; it does not establish that the parameter is unimportant in every setting.

A large $\sigma_i$ means effects vary across sampled locations. This can arise from nonlinearity, interactions with other inputs, or a change in the sign of the effect. It can also arise from an inadequately sampled trajectory set or numerical noise. Therefore, $\sigma_i$ is a diagnostic clue, not a unique measurement of interaction or nonlinearity. A large $\mu_i^*$ together with a large $\sigma_i$ identifies a parameter whose effects may be substantial but whose response is not well described by one consistent local slope.

Worked Example: Summarizing Effects for One PBPK Parameter

Suppose a hypothetical screening run uses four trajectories for a parameter that affects predicted AUC. The calculated elementary effects are 2, 3, $-1$, and 4 mg·h/L per normalized input unit. Calculate the three summaries.

The signed mean is $$ \mu_i=\frac{2+3-1+4}{4}=\frac{8}{4}=2\text{ mg·h/L per normalized input unit}. $$

The mean absolute effect is $$ \mu_i^*=\frac{|2|+|3|+|-1|+|4|}{4} =\frac{10}{4}=2.5\text{ mg·h/L per normalized input unit}. $$

For the sample standard deviation, deviations from the signed mean of 2 are $0,1,-3,2$. Thus $$ \sigma_i=\sqrt{\frac{0^2+1^2+(-3)^2+2^2}{4-1}} =\sqrt{\frac{14}{3}} =2.160\text{ mg·h/L per normalized input unit}. $$

The positive $\mu_i$ indicates a positive average signed effect, but one trajectory has a negative effect. The larger $\mu_i^*$ and the substantial $\sigma_i$ retain information about the magnitude and variability that the signed mean alone would obscure. With only four trajectories, these are preliminary screening summaries, not precise estimates of every possible effect in the parameter space.

Why Effects Vary: A Short Derivation

The connection between elementary effects and local derivatives follows from a Taylor expansion. For a fixed setting of all other inputs, suppose $f$ has a bounded third derivative with respect to $z_i$ in a neighborhood of the point. Then a forward step gives

f(\mathbf{z}+\Delta e_i) =f(\mathbf{z})+\Delta\frac{\partial f}{\partial z_i}(\mathbf{z}) +\frac{\Delta^2}{2}\frac{\partial^2 f}{\partial z_i^2}(\mathbf{z}) +O(\Delta^3).

Subtract $f(\mathbf{z})$ and divide by $\Delta$. Since $\Delta$ is nonzero, the result is

EE_i(\mathbf{z})= \frac{\partial f}{\partial z_i}(\mathbf{z}) +\frac{\Delta}{2}\frac{\partial^2 f}{\partial z_i^2}(\mathbf{z}) +O(\Delta^2).

Thus an elementary effect approximates a local derivative, with a finite-step contribution from curvature. The Morris design samples such effects at multiple settings of the other inputs. If the slope changes with those settings, the effects vary; curvature in $z_i$ can also make the finite-step effect differ from the derivative at the starting point. This explains why $\sigma_i$ can flag complex response behavior, but cannot by itself distinguish its cause. The stated remainder requires the stronger bounded-third-derivative condition; twice differentiability alone would not justify an $O(\Delta^3)$ remainder in the first expansion.

Worked Example: An Interaction Makes Effects Depend on Location

Consider a hypothetical normalized response $$ f(z_1,z_2)=2z_1+3z_1z_2, $$ where $z_1$ and $z_2$ could represent normalized PBPK parameters. With $p=4$, use $\Delta=2/3$. Suppose trajectories sample the effect for $z_1$ at $z_2=0$ and $z_2=1$, with valid starting values for $z_1$.

At $z_2=0$, the effect for $z_1$ is $$ \frac{f(z_1+2/3,0)-f(z_1,0)}{2/3} =\frac{2(z_1+2/3)-2z_1}{2/3}=2. $$ At $z_2=1$, it is $$ \frac{f(z_1+2/3,1)-f(z_1,1)}{2/3} =\frac{5(z_1+2/3)-5z_1}{2/3}=5. $$ Therefore, for these two effects, $\mu_1=\mu_1^*=(2+5)/2=3.5$, and $$ \sigma_1=\sqrt{\frac{(2-3.5)^2+(5-3.5)^2}{2-1}} =\sqrt{4.5}=2.121. $$

For $z_2$, sample effects at $z_1=0$ and $z_1=1$. The corresponding effects are 0 and 3, since changing $z_2$ changes only the interaction term $3z_1z_2$. Hence $\mu_2=\mu_2^*=1.5$ and $\sigma_2=\sqrt{[(0-1.5)^2+(3-1.5)^2]/1}=2.121$. The interaction makes each parameter's effect depend on the value of the other. This example is linear in each coordinate when the other is fixed, so the variation here is due to interaction, not curvature in an individual coordinate.

Using Results to Guide Screening

A useful first display is a plot of $\mu_i^*$ against $\sigma_i$, with one point per parameter and units and endpoint stated. Parameters with small $\mu_i^*$ over the chosen ranges are candidates for lower priority in follow-up analysis, whereas those with large $\mu_i^*$ merit attention. A large $\sigma_i$ cautions against describing a parameter by a single constant effect. This display supports prioritization, but does not alone justify fixing or removing a parameter from a calibrated PBPK model.

Worked Example: Comparing Two Hypothetical Parameters

For a hypothetical AUC endpoint in mg·h/L, four trajectories give effects of $0.12, 0.10, 0.14,$ and $0.11$ for parameter A. Parameter B has effects of $-0.8, 0.7, -0.9,$ and $0.8$. Calculate the mean absolute effects and sample standard deviations.

For A, $$ \mu_A^*=\frac{0.12+0.10+0.14+0.11}{4}=0.1175\text{ mg·h/L per normalized input unit}. $$ Its signed mean is $0.1175$. The squared deviations from that mean sum to $0.000875$, so $$ \sigma_A=\sqrt{\frac{0.000875}{3}}=0.01708\text{ mg·h/L per normalized input unit}. $$

For B, $$ \mu_B^*=\frac{0.8+0.7+0.9+0.8}{4}=0.8\text{ mg·h/L per normalized input unit}, $$ while its signed mean is $(-0.8+0.7-0.9+0.8)/4=-0.05$. The squared deviations from $-0.05$ sum to $2.57$, giving $$ \sigma_B=\sqrt{\frac{2.57}{3}}=0.9256\text{ mg·h/L per normalized input unit}. $$

In this hypothetical screen, B has a much larger average effect magnitude and much greater effect variability. Its near-zero signed mean is misleading because positive and negative effects largely cancel. The results suggest investigating B across the parameter space; they do not establish which mechanism causes the variation.

Assumptions, Diagnostics, and Alternatives

Morris screening is most straightforward when the model output is a reproducible scalar function of the inputs, the selected parameter ranges are meaningful, and the sampled combinations are physiologically interpretable. Standard trajectories change one normalized coordinate at a time while holding the others fixed. This construction is conceptually clean when parameters can be varied independently over the defined ranges. If parameter dependence or physiological constraints make many combinations implausible, unconstrained trajectories can evaluate irrelevant model states; the resulting summaries then need careful qualification.

Check that each simulation completes and that the output is on the intended scale. Repeat the screen with more trajectories and compare rankings and effect distributions. If rankings shift substantially, the initial trajectory sample may be too small. Inspect effects for signs, extreme values, and clusters rather than relying only on the three summaries. When an effect is small relative to known numerical solver error or simulation variability, it may not represent a reliable parameter response. Tightening numerical tolerances or assessing replicate simulation noise can help separate numerical instability from genuine effect variation.

Morris screening differs from local sensitivity analysis, which evaluates derivatives or finite differences near a specified parameter vector. The local approach can be more directly tied to that reference condition, but may not reflect behavior across broad ranges. Morris screening samples multiple locations and is often more economical than a dense multidimensional grid, but its summaries are approximate and depend on the selected trajectories. Neither method removes the need to choose relevant endpoints and ranges. Morris summaries also do not allocate a fraction of total output variance to each parameter; they are elementary-effect summaries, not a variance decomposition.

A common mistake is to interpret a large $\sigma_i$ as proof of a specific interaction. In fact, high variability may reflect interactions, nonlinear response, sign changes, sparse trajectory coverage, or numerical noise. Another mistake is to interpret a small $\mu_i^*$ as proof that the parameter can safely be fixed. The correction is to treat the result as conditional screening evidence: state the output and ranges, inspect diagnostics, and verify any consequential modeling decision with additional analysis.

Key takeaway: Morris screening uses a modest number of structured trajectories to summarize how finite parameter changes affect a chosen PBPK output across sampled locations. Use $\mu_i^*$ to describe average effect magnitude and $\sigma_i$ to flag variability, while treating both as range- and design-dependent screening evidence.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions.

A parameter has effects $-2$, 2, and $-2$ in three trajectories. What are its signed mean and mean absolute effect? What does their difference suggest?
For a PBPK endpoint measured in mg/L, what are the units of an elementary effect when its parameter input is normalized to the interval from 0 to 1?
Why is a high $\sigma_i$ not sufficient evidence to conclude that parameter $i$ interacts with another parameter?
Give two reasons a parameter with a small $\mu_i^*$ in one screen might still warrant further examination.

Solutions. (1) The signed mean is $(-2+2-2)/3=-2/3$; the mean absolute effect is $(2+2+2)/3=2$. The difference indicates cancellation in the signed mean and effects that change sign. (2) The effect has units of mg/L per normalized input unit. (3) Variability can also arise from curvature, sign changes, limited sampling, or numerical noise. (4) The chosen ranges or endpoint may not represent the intended application, and the parameter could be important under another dose, scenario, output, or constrained combination of inputs.