Worked Screening Analysis for a Cytokine Pathway Model

Screening a Cytokine Pathway Model

A mechanistic model may contain many parameters, while only a few strongly influence an output that matters for a particular question. Morris screening is a computationally economical way to identify parameters whose changes produce substantial output changes over specified ranges. It can also flag parameters whose effects vary across the sampled parameter space, a pattern that may arise from nonlinearity or interactions.

This tutorial works through a hypothetical cytokine production-and-clearance model. The target is the steady-state cytokine concentration under a fixed stimulus. This is a deliberate output choice: it is relevant when the scientific question concerns sustained stimulation, and it yields a transparent worked calculation. A transient response, peak concentration, or threshold-crossing time would be a different target and could produce a different screening order.

Model, Target, and Screening Assumptions

Let \(C(t)\) be cytokine concentration in ng/mL. Under a constant dimensionless stimulus \(S\), suppose its dynamics are

Here \(q_{\max}\) is the maximum production rate in ng/mL/h, \(EC_{50}\) is the stimulus level giving half-maximal production, and \(k_{\mathrm{clear}}\) is the first-order clearance rate in h\(^{-1}\). For this worked analysis, fix \(S=1\), assume \(k_{\mathrm{clear}}>0\), and use the equilibrium concentration as the scalar output.

The screening calculation relies on a deterministic model, defined parameter ranges, a fixed output definition, and parameter values for which the model is well-defined. The trajectory design must stay within the ranges. Here, parameters are varied one at a time along each trajectory; each trajectory uses the same step size in normalized coordinates. Morris summaries describe effects over the selected ranges and design points, not probabilities that parameters are important.

The equilibrium target follows directly from the model. Setting \(dC/dt=0\) and using \(S=1\) gives \(q_{\max}/(EC_{50}+1)=k_{\mathrm{clear}}C^*\), so

This is a stable equilibrium: for a constant production rate, the difference from equilibrium decays at rate \(k_{\mathrm{clear}}\). In fact, \(C(t)=C^*+[C(0)-C^*]e^{-k_{\mathrm{clear}}t}\). Thus the target represents sustained exposure to the fixed stimulus, not an early transient response.

Ranges and Trajectory Design

Use the following hypothetical ranges. They are chosen to make the calculations explicit, not asserted as empirical ranges for a particular cytokine.

Each normalized coordinate runs from 0 to 1. With three factors and Morris step \(\Delta=2/3\), grid points are 0, 1/3, 2/3, and 1; a valid step joins grid points separated by \(2/3\). For example, a step in \(u_q\) from 0 to \(2/3\) changes \(q_{\max}\) from 1 to \(1+2(2/3)=7/3\) ng/mL/h. An effect is therefore expressed in ng/mL per unit change in the normalized factor. To report change per physical unit of \(q_{\max}\), divide the normalized-coordinate effect by the \(q_{\max}\) range width of 2 ng/mL/h.

A Morris trajectory changes each factor once. The two trajectories below are hypothetical and illustrate calculations; two paths alone are not enough to establish a stable ranking for a real analysis. For each path, the factor order is \(q_{\max}\), \(k_{\mathrm{clear}}\), then \(EC_{50}\).

Summarizing Effects and Interpreting Them

For \(r\) trajectories, a common Morris summary for factor \(i\) is the mean absolute elementary effect, \(\mu_i^*\), together with the sample standard deviation \(s_i\) of its signed elementary effects. The mean absolute effect avoids cancellation between positive and negative effects; the standard deviation describes how much signed effects vary across the sampled trajectories.

A useful justification for retaining the sign in the denominator follows from a local expansion. If the output is differentiable in normalized coordinate \(u_i\), then for a small signed step \(h\),

The signed quotient therefore approximates the derivative in the direction of increasing \(u_i\), whether the actual path step is positive or negative. Dividing by \(|h|\) instead would reverse the interpretation on a decreasing step. Also, a large \(s_i\) means effects differ across sampled locations; it can reflect curvature, interactions, or both, and by itself does not identify which explanation is responsible.

Diagnostics, Alternatives, and Follow-Up

Before interpreting a ranking, verify that the model and design behaved as intended. In this example, check that every parameter lies within its stated range, \(k_{\mathrm{clear}}\) remains positive, the numerical solver reaches the intended equilibrium if simulation rather than the closed-form expression is used, and the output has consistent units. For more complex models, inspect solver warnings, failed runs, implausible state values, and whether small tolerances or alternative initial conditions materially change the selected output.

A simple analytic check is available here. The equilibrium increases with \(q_{\max}\), decreases with \(k_{\mathrm{clear}}\), and decreases with \(EC_{50}\); the signs of the computed effects agree. A sign disagreement can be a warning of a calculation, coding, or parameter-mapping error, although nonlinear models can have effects that change sign over their ranges. Plotting output against each parameter at several settings of the other parameters can help reveal that pattern.

Morris screening is a qualitative-to-semiquantitative screen, not a replacement for every global sensitivity method. Local derivatives are cheaper and can be informative near one calibrated parameter vector, but do not represent variation across broad ranges. Variance-based indices such as Sobol indices quantify contributions to output variance under specified input distributions, but generally require many more model evaluations. Morris elementary effects instead use selected grid trajectories and can expose broad-range influence at lower cost. The methods answer different questions, so a large Morris effect is not itself a variance fraction or a causal conclusion about biological intervention.

Common failure modes include choosing ranges that include invalid parameter values, screening an output that does not match the scientific objective, and treating a small number of paths as conclusive. For time-course models, selecting only one time point can miss a parameter that matters strongly at another time or changes the shape of the response. It is also a mistake to average signed effects and infer low importance when positive and negative effects cancel; use \(\mu_i^*\) for magnitude and inspect signed effects to understand direction. Conversely, large \(s_i\) should prompt investigation, not an automatic claim of interaction.