---
title: "Effect Compartment Pharmacodynamic Models"
execute:
  cache: true
  error: false
engine: julia
author:
  - Vijay Ivaturi
format:
  html:
    include-in-header:
      text: |
        <style>
          .header-container {
            width: 100%;
          }
          .logo {
            display: block;
            width: 150px;
            margin-left: -2.5%;
          }
        </style>
        <link rel = "shortcut icon" href = "https://pumas.ai/favicon.ico" />
    include-before-body:
      text: |
        <div class="header-container">
          <img src="https://pumas-assets.s3.amazonaws.com/CompanyLogos/PumasAI/RGB/PNG/Pumas+AI+Primary%404x.png" alt="Pumas Logo" class="logo">
        </div>
    self-contained-math: true
    html-math-method: katex
    anchor-sections: true
    theme: default
    toc: true
    toc-depth: 4
    toc-expand: 2
    toc-location: left
    toc-title: Contents
    number-sections: true
    code-summary: Show/Hide Code
fig-format: svg
fig-width: 8
fig-height: 6
license: CC BY-SA 4.0
bibliography: references.bib
mermaid:
  theme: neutral
---

```{julia}
#| echo: false
#| output: false
using Pumas
using PumasUtilities
using PumasPlots
using CairoMakie
using AlgebraOfGraphics
using DataFrames
using DataFramesMeta
using Random
```

## Learning Objectives

By the end of this tutorial, you will be able to:

  - Recognize hysteresis in concentration-effect data
  - Understand the effect compartment (biophase) concept
  - Implement effect compartment models in Pumas
  - Interpret the $k_{e0}$ parameter and equilibration half-life
  - Compare effect compartment models with direct response models
  - Identify appropriate use cases for this model type

## Introduction

When plotting drug effect against plasma concentration,
we sometimes observe that the curve forms a **loop** rather than overlapping.
This phenomenon is called **hysteresis** and indicates a temporal disconnect
between concentration and effect.

```{mermaid}
%%| label: fig-effect-compartment-concept
%%| fig-cap: "Effect compartment model: Effect site concentration (Ce) lags behind plasma concentration (Cp)"
flowchart LR
    subgraph PK["Pharmacokinetics"]
        Dose[Drug Dose] --> Cp[Plasma<br/>Concentration Cp]
    end
    subgraph Biophase["Effect Site (Biophase)"]
        Cp -->|ke0| Ce[Effect Site<br/>Concentration Ce]
    end
    subgraph PD["Pharmacodynamics"]
        Ce --> Effect["Effect<br/>E = f(Ce)"]
    end

    style Ce fill:#fff3e0
```

### When Direct Response Fails

Consider a patient receiving morphine for pain relief.
You measure plasma concentrations and pain scores over time:

  - At 1 hour: Plasma concentration is 30 ng/mL, pain score is 7
  - At 3 hours: Plasma concentration is 20 ng/mL, pain score is 3
  - At 6 hours: Plasma concentration is 10 ng/mL, pain score is 5

If you plot pain relief vs. concentration, the points don't follow a single line -
they form a loop. This is counter-clockwise hysteresis:
the effect is lower at the same concentration during the rising phase
compared to the falling phase.

:::{.callout-note}

## Types of Hysteresis

**Counter-clockwise (proteresis)**: Effect lags behind concentration.
The drug takes time to distribute to the effect site.
Modeled with an effect compartment.

**Clockwise**: Effect precedes or exceeds what concentration would predict.
Usually indicates tolerance, receptor down-regulation, or indirect mechanisms.
Typically requires different modeling approaches.
:::

## The Effect Compartment Concept

@sheiner1979 introduced the **effect compartment** (also called **biophase**)
to explain temporal delays between plasma concentration and drug effect.

### Key Principles

 1. The effect site is a hypothetical compartment linked to the central (plasma) compartment
 2. Drug equilibrates between plasma and effect site with first-order kinetics
 3. The effect compartment is assumed to have negligible volume (doesn't affect PK)
 4. Effect is driven by the **effect site concentration** ($C_e$), not plasma concentration ($C_p$)

### Mathematical Framework

The effect compartment is described by a first-order differential equation:

$$
\frac{dC_e}{dt} = k_{e0} \cdot (C_p - C_e)
$$

Where:

| Parameter | Description                             | Units       |
|:--------- |:--------------------------------------- |:----------- |
| $C_e$     | Effect site concentration               | mass/volume |
| $C_p$     | Plasma concentration                    | mass/volume |
| $k_{e0}$  | First-order equilibration rate constant | 1/time      |

### Equilibration Half-Life

The $k_{e0}$ parameter has a direct clinical interpretation:

$$
t_{1/2,e0} = \frac{\ln(2)}{k_{e0}} = \frac{0.693}{k_{e0}}
$$

This is the **equilibration half-life** - the time for the effect site concentration
to reach 50% of the plasma concentration after a step change.

| $k_{e0}$ (1/hr) | $t_{1/2,e0}$ (hr) | Interpretation         |
|:--------------- |:----------------- |:---------------------- |
| 0.693           | 1.0               | Fast equilibration     |
| 0.231           | 3.0               | Moderate equilibration |
| 0.0693          | 10.0              | Slow equilibration     |

:::{.callout-tip}

## Clinical Interpretation

A drug with $t_{1/2,e0}$ = 30 minutes will reach ~94% equilibration
between plasma and effect site within 2 hours (4 half-lives).
:::

## Pumas Implementation

### Effect Compartment Model

```{julia}
effect_compartment_model = @model begin
    @metadata begin
        desc = "One-Compartment PK with Effect Compartment PD Model"
        timeu = u"hr"
    end

    @param begin
        # PK parameters
        """
        Absorption rate constant (1/hr)
        """
        pop_Ka ∈ RealDomain(; lower = 0, init = 1.0)
        """
        Clearance (L/hr)
        """
        pop_CL ∈ RealDomain(; lower = 0, init = 0.5)
        """
        Volume of distribution (L)
        """
        pop_Vc ∈ RealDomain(; lower = 0, init = 20.0)

        # PD parameters
        """
        Baseline response
        """
        pop_E0 ∈ RealDomain(; lower = 0, init = 30.0)
        """
        Maximum drug-induced effect
        """
        pop_Emax ∈ RealDomain(; lower = 0, init = 80.0)
        """
        Effect site concentration at 50% Emax
        """
        pop_EC50 ∈ RealDomain(; lower = 0, init = 1.0)
        """
        Effect site equilibration rate constant (1/hr)
        """
        pop_ke0 ∈ RealDomain(; lower = 0, init = 0.5)

        # Between-subject variability
        Ω_pk ∈ PDiagDomain(3)
        Ω_pd ∈ PDiagDomain(2)

        # Residual unexplained variability
        σ_prop_pk ∈ RealDomain(; lower = 0, init = 0.1)
        σ_add_pd ∈ RealDomain(; lower = 0, init = 3.0)
    end

    @random begin
        η_pk ~ MvNormal(Ω_pk)
        η_pd ~ MvNormal(Ω_pd)
    end

    @pre begin
        # Individual PK parameters
        CL = pop_CL * exp(η_pk[1])
        Vc = pop_Vc * exp(η_pk[2])
        Ka = pop_Ka * exp(η_pk[3])

        # Individual PD parameters
        E0 = pop_E0 * exp(η_pd[1])
        Emax = pop_Emax * exp(η_pd[2])
        EC50 = pop_EC50
        ke0 = pop_ke0
    end

    @init begin
        # Effect compartment starts at zero (no drug initially)
        Ce = 0.0
    end

    @vars begin
        # Effect site concentration for output
        Ceff = Ce
        # Plasma concentration (for plotting)
        Cp = Central / Vc
        # Predicted effect (for plotting)
        Effect = E0 + (Emax * Ce) / (EC50 + Ce)
    end

    @dynamics begin
        # PK dynamics
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central

        # Effect compartment dynamics
        # Ce equilibrates toward Cp with rate ke0
        Ce' = ke0 * (Central / Vc - Ce)
    end

    @derived begin
        # Observations with residual error
        conc ~ @. Normal(Cp, abs(Cp) * σ_prop_pk)
        effect ~ @. Normal(Effect, σ_add_pd)
    end
end
```

### Understanding the Model Structure

:::{.panel-tabset}

#### Effect Compartment Dynamics

The effect compartment dynamics are key:

```julia
Ce' = ke0 * (Central / Vc - Ce)
```

This says:

  - $C_e$ changes proportionally to the difference between $C_p$ and $C_e$
  - When $C_p > C_e$: $C_e$ increases (drug distributing into effect site)
  - When $C_p < C_e$: $C_e$ decreases (drug redistributing out)
  - Rate is governed by $k_{e0}$

#### Initial Conditions

The `@init` block sets $C_e = 0$ at time zero.
This assumes no drug at the effect site before dosing.

#### Effect Calculation

The effect is calculated from $C_e$, not $C_p$:

```julia
Effect := @. E0 + (Emax * Ce) / (EC50 + Ce)
```

This is why the effect lags behind concentration.

:::

## Simulation Workflow

### Comparing Cp and Ce

Let's visualize how the effect site concentration lags behind plasma concentration:

```{julia}
# Parameters with moderate ke0 (t1/2,e0 ≈ 1.4 hr)
params_ec = (
    pop_Ka = 1.5,
    pop_CL = 0.5,
    pop_Vc = 20.0,
    pop_E0 = 30.0,
    pop_Emax = 80.0,
    pop_EC50 = 1.0,
    pop_ke0 = 0.5,  # t1/2,e0 = 0.693/0.5 ≈ 1.4 hr
    Ω_pk = Diagonal([0.04, 0.04, 0.04]),
    Ω_pd = Diagonal([0.04, 0.04]),
    σ_prop_pk = 0.1,
    σ_add_pd = 3.0,
)

# Dosing regimen
dosing = DosageRegimen(100, time = 0, cmt = 1)

# Single subject
subject = Subject(id = 1, events = dosing)

# Simulate (without residual error for clean visualization)
sim_ec = simobs(
    effect_compartment_model,
    subject,
    params_ec,
    obstimes = 0:0.25:24,
    simulate_error = false,
)
```

```{julia}
#| label: fig-cp-ce-comparison
#| fig-cap: "Comparison of plasma concentration (Cp) and effect site concentration (Ce)"

# Convert simulation to DataFrame
sim_df = DataFrame(sim_ec)
dropmissing!(sim_df, [:Cp, :Ceff])  # Ensure no missing values for plotting
# Reshape data for plotting Cp and Ce together
conc_df = stack(
    select(sim_df, :time, :Cp, :Ceff),
    [:Cp, :Ceff];
    variable_name = :concentration_type,
    value_name = :concentration,
)
conc_df.concentration_type =
    replace.(
        string.(conc_df.concentration_type),
        "Cp" => "Cp (plasma)",
        "Ceff" => "Ce (effect site)",
    )

# Create plot
plt =
    data(conc_df) *
    mapping(
        :time => "Time (hr)",
        :concentration => "Concentration (mg/L)",
        color = :concentration_type => "Type",
    ) *
    visual(Lines, linewidth = 2)

# Find Tmax for each
cp_max_idx = argmax(sim_df.Cp)
ce_max_idx = argmax(sim_df.Ceff)
tmax_df = DataFrame(
    time = [sim_df.time[cp_max_idx], sim_df.time[ce_max_idx]],
    concentration = [sim_df.Cp[cp_max_idx], sim_df.Ceff[ce_max_idx]],
    concentration_type = ["Cp (plasma)", "Ce (effect site)"],
)

tmax_plt =
    data(tmax_df) *
    mapping(:time, :concentration, color = :concentration_type => "Type") *
    visual(Scatter, markersize = 12)

draw(
    plt + tmax_plt;
    axis = (title = "Plasma vs Effect Site Concentration", width = 700, height = 400),
)
```

:::{.callout-note title="Key Observation"}
Notice that $C_e$ peaks later than $C_p$, is always lower than $C_p$ during the absorption phase, and can exceed $C_p$ during the elimination phase (when $C_p$ is falling faster).
:::

### Visualizing Hysteresis

The hallmark of an effect compartment model is hysteresis in the effect-concentration plot:

```{julia}
#| label: fig-hysteresis
#| fig-cap: "Hysteresis loop: Effect vs concentration shows counter-clockwise pattern"

# Use sim_df from previous code block (Cp and Effect are in @vars, no missing values)
hysteresis_df = select(sim_df, :Cp, :Effect)
conc = hysteresis_df.Cp
eff = hysteresis_df.Effect

# Create base hysteresis plot
plt =
    data(hysteresis_df) *
    mapping(:Cp => "Plasma Concentration (mg/L)", :Effect => "Effect") *
    visual(Lines, linewidth = 2, color = :purple)

fig = Figure(size = (600, 500))
ax = Axis(
    fig[1, 1],
    xlabel = "Plasma Concentration (mg/L)",
    ylabel = "Effect",
    title = "Effect vs Plasma Concentration (Hysteresis)",
)
draw!(ax, plt)

# Add arrows to show direction (every 4th point)
for i = 4:4:length(conc)-1
    dx = conc[i+1] - conc[i]
    dy = eff[i+1] - eff[i]
    arrows!(
        ax,
        [conc[i]],
        [eff[i]],
        [dx * 0.5],
        [dy * 0.5],
        color = :purple,
        arrowsize = 10,
    )
end

# Mark start and end
scatter!(
    ax,
    [conc[1]],
    [eff[1]],
    markersize = 15,
    color = :green,
    marker = :circle,
    label = "Start (t=0)",
)
scatter!(
    ax,
    [conc[end]],
    [eff[end]],
    markersize = 15,
    color = :red,
    marker = :star5,
    label = "End (t=24h)",
)

# Add text for loop direction
text!(ax, 2.5, 75, text = "Counter-clockwise\n(proteresis)", fontsize = 14)

axislegend(ax, position = :rb)
fig
```

### Impact of ke0 on Delay

Let's compare different $k_{e0}$ values to see how they affect the delay:

```{julia}
#| label: fig-ke0-comparison
#| fig-cap: "Effect of ke0 on the temporal relationship between Cp and effect"

ke0_values = [0.2, 0.5, 2.0]
labels = ["ke0 = 0.2 (slow)", "ke0 = 0.5 (moderate)", "ke0 = 2.0 (fast)"]

# Simulate for each ke0 value and collect results (no residual error for visualization)
ke0_results = map(zip(ke0_values, labels)) do (ke0, label)
    params_var = merge(params_ec, (; pop_ke0 = ke0))
    sim_var = simobs(
        effect_compartment_model,
        subject,
        params_var,
        obstimes = 0:0.25:24,
        simulate_error = false,
    )
    df = DataFrame(sim_var)
    df.ke0_label .= label
    df
end
ke0_df = vcat(ke0_results...)

# Time course plot
time_plt =
    data(ke0_df) *
    mapping(:time => "Time (hr)", :Effect => "Effect", color = :ke0_label => "ke0") *
    visual(Lines, linewidth = 2)

# Hysteresis plot
hyst_plt =
    data(ke0_df) *
    mapping(
        :Cp => "Plasma Concentration",
        :Effect => "Effect",
        color = :ke0_label => "ke0",
    ) *
    visual(Lines, linewidth = 2)

fig = Figure(size = (900, 400))
draw!(fig[1, 1], time_plt; axis = (title = "Effect Time Course for Different ke0",))
draw!(fig[1, 2], hyst_plt; axis = (title = "Hysteresis Loops for Different ke0",))
fig
```

:::{.callout-tip title="Interpretation"}
**Low $k_{e0}$** (0.2): Large hysteresis loop, effect lags significantly.
**High $k_{e0}$** (2.0): Small loop, approaches direct response behavior.
As $k_{e0} \to \infty$: No hysteresis, $C_e$ = $C_p$ (direct response).
:::

## Comparison with Direct Response

Let's directly compare the effect compartment model with a direct response model:

```{julia}
# Direct response model (ke0 → infinity, effectively Cp = Ce)
direct_model = @model begin
    @metadata begin
        desc = "Direct Response for Comparison"
        timeu = u"hr"
    end

    @param begin
        pop_Ka ∈ RealDomain(; lower = 0)
        pop_CL ∈ RealDomain(; lower = 0)
        pop_Vc ∈ RealDomain(; lower = 0)
        pop_E0 ∈ RealDomain(; lower = 0)
        pop_Emax ∈ RealDomain(; lower = 0)
        pop_EC50 ∈ RealDomain(; lower = 0)
        Ω_pk ∈ PDiagDomain(3)
        Ω_pd ∈ PDiagDomain(2)
        σ_prop_pk ∈ RealDomain(; lower = 0)
        σ_add_pd ∈ RealDomain(; lower = 0)
    end

    @random begin
        η_pk ~ MvNormal(Ω_pk)
        η_pd ~ MvNormal(Ω_pd)
    end

    @pre begin
        CL = pop_CL * exp(η_pk[1])
        Vc = pop_Vc * exp(η_pk[2])
        Ka = pop_Ka * exp(η_pk[3])
        E0 = pop_E0 * exp(η_pd[1])
        Emax = pop_Emax * exp(η_pd[2])
        EC50 = pop_EC50
    end

    @dynamics begin
        Depot' = -Ka * Depot
        Central' = Ka * Depot - (CL / Vc) * Central
    end

    @vars begin
        Cp = Central / Vc
        Effect = E0 + (Emax * Cp) / (EC50 + Cp)
    end

    @derived begin
        conc ~ @. Normal(Cp, abs(Cp) * σ_prop_pk)
        effect ~ @. Normal(Effect, σ_add_pd)
    end
end

params_direct = (
    pop_Ka = 1.5,
    pop_CL = 0.5,
    pop_Vc = 20.0,
    pop_E0 = 30.0,
    pop_Emax = 80.0,
    pop_EC50 = 1.0,
    Ω_pk = Diagonal([0.04, 0.04, 0.04]),
    Ω_pd = Diagonal([0.04, 0.04]),
    σ_prop_pk = 0.1,
    σ_add_pd = 3.0,
)

sim_direct = simobs(
    direct_model,
    subject,
    params_direct,
    obstimes = 0:0.25:24,
    simulate_error = false,
)
```

```{julia}
#| label: fig-direct-vs-ec
#| fig-cap: "Comparison of direct response and effect compartment models"

# Convert both simulations to DataFrames
direct_df = DataFrame(sim_direct)
direct_df.model .= "Direct Response"
ec_df = DataFrame(sim_ec)
ec_df.model .= "Effect Compartment"

# Combine for comparison (using vars: Cp and Effect)
comparison_df = vcat(
    select(direct_df, :time, :Cp, :Effect, :model),
    select(ec_df, :time, :Cp, :Effect, :model),
)

# Time course plot
time_plt =
    data(comparison_df) *
    mapping(:time => "Time (hr)", :Effect => "Effect", color = :model => "Model") *
    visual(Lines, linewidth = 2)

# Concentration-effect plot
ce_plt =
    data(comparison_df) *
    mapping(:Cp => "Plasma Concentration", :Effect => "Effect", color = :model => "Model") *
    visual(Lines, linewidth = 2)

fig = Figure(size = (900, 400))
draw!(fig[1, 1], time_plt; axis = (title = "Effect Time Course",))
draw!(fig[1, 2], ce_plt; axis = (title = "Concentration-Effect Relationship",))
fig
```

## Population Simulation

```{julia}
# Create population
population = map(1:20) do i
    Subject(id = i, events = dosing)
end

# Simulate
Random.seed!(123)
sim_pop = simobs(effect_compartment_model, population, params_ec, obstimes = 0:0.5:24)
```

```{julia}
#| label: fig-population-ec
#| fig-cap: "Population simulation showing variability in effect compartment model"

# Convert population simulation to DataFrame
pop_df = DataFrame(sim_pop)

# Time course plot (using vars: Cp, Effect - always available for each time point)
time_pop_plt =
    data(pop_df) *
    mapping(:time => "Time (hr)", :Effect => "Effect", group = :id) *
    visual(Lines, linewidth = 1, alpha = 0.5, color = :coral)

# Hysteresis plot
hyst_pop_plt =
    data(pop_df) *
    mapping(:Cp => "Plasma Concentration", :Effect => "Effect", group = :id) *
    visual(Lines, linewidth = 1, alpha = 0.5, color = :coral)

fig = Figure(size = (900, 400))
draw!(fig[1, 1], time_pop_plt; axis = (title = "Population Effect Profiles",))
draw!(fig[1, 2], hyst_pop_plt; axis = (title = "Population Hysteresis",))
fig
```

## Fitting Workflow

### Generate Synthetic Data

```{julia}
#| output: false
Random.seed!(456)
sim_data = simobs(
    effect_compartment_model,
    population,
    params_ec,
    obstimes = [0.5, 1, 2, 4, 6, 8, 12, 24],
)
df_sim = DataFrame(sim_data)
```

### Fit the Model

```{julia}
fit_pop = read_pumas(
    df_sim;
    id = :id,
    time = :time,
    amt = :amt,
    evid = :evid,
    cmt = :cmt,
    observations = [:conc, :effect],
)

init_params = (
    pop_Ka = 1.2,
    pop_CL = 0.4,
    pop_Vc = 25.0,
    pop_E0 = 25.0,
    pop_Emax = 70.0,
    pop_EC50 = 0.8,
    pop_ke0 = 0.4,
    Ω_pk = Diagonal([0.02, 0.02, 0.02]),
    Ω_pd = Diagonal([0.02, 0.02]),
    σ_prop_pk = 0.15,
    σ_add_pd = 5.0,
)

fit_result = fit(effect_compartment_model, fit_pop, init_params, FOCE())
```

```{julia}
# Estimated parameters
coef(fit_result)
```

## Model Diagnostics

```{julia}
insp = inspect(fit_result)
```

```{julia}
#| label: fig-gof-ec
#| fig-cap: "Goodness-of-fit for effect compartment model"
goodness_of_fit(insp; observations = [:effect])
```

## Clinical Examples

:::{.panel-tabset}

### Morphine

  - **Drug**: Morphine (analgesic)
  - **Site of action**: Central nervous system
  - **$t_{1/2,e0}$**: ~1.5-2 hours
  - **Clinical implication**: Pain relief lags behind plasma levels; titration requires patience
  - **Reference**: @lotsch2006

### Digoxin

  - **Drug**: Digoxin (cardiac glycoside)
  - **Site of action**: Cardiac muscle
  - **$t_{1/2,e0}$**: ~8 hours
  - **Clinical implication**: Slow equilibration means toxicity may appear delayed
  - **Reference**: @jelliffe2014

### Warfarin

  - **Drug**: Warfarin (anticoagulant)
  - **Site of action**: Liver (clotting factor synthesis)
  - **$t_{1/2,e0}$**: ~6-8 hours for initial effects
  - **Note**: Full anticoagulant effect involves indirect response (clotting factor turnover)
  - **Reference**: @hamberg2007

:::

## Summary

In this tutorial, we covered:

 1. **Hysteresis**: Counter-clockwise loops in effect-concentration plots indicate effect site distribution delay
 2. **Effect compartment**: A hypothetical compartment linked to plasma that explains the delay
 3. **$k_{e0}$ parameter**: The equilibration rate constant; $t_{1/2,e0} = 0.693/k_{e0}$
 4. **Pumas implementation**: Using the `@dynamics` block for $C_e$ dynamics
 5. **Comparison**: Effect compartment → direct response as $k_{e0}$ increases

### When to Use Effect Compartment Models

  - Counter-clockwise hysteresis is observed
  - Peak effect occurs after peak concentration
  - Drug distributes slowly to site of action
  - CNS-acting drugs with blood-brain barrier

### Next Steps

If your data shows:

  - **Effect persists long after drug clears**: Consider [Indirect Response Models](03-indirect-response.qmd)
  - **Irreversible drug-target interaction**: Consider [K-PD Models](04-kpd-models.qmd)
  - **Model selection challenges**: See [Model Selection Tutorial](05-model-selection.qmd)
