T5 - Using interventions

In this tutorial we build a screen-and-treat program for Nigeria — screening, triage, and two treatment options wired together — and measure how many cervical cancers it prevents at two levels of coverage. We then add therapeutic vaccination on top of it.

Products and programs

HPVsim separates the product from the program that delivers it.

A product is a test, a treatment or a vaccine, with its own performance. The ones that ship with HPVsim are in three CSV files you can read directly:

screening and triage tests — VIA, cytology, HPV DNA, HPV16/18 typing — in hpvsim/data/products_dx.csv
treatments — ablation, excision, radiation, and the therapeutic vaccines — in hpvsim/data/products_tx.csv
prophylactic vaccines — bivalent, quadrivalent, nonavalent — in hpvsim/data/products_vx.csv

A program decides who is offered the product, when, and how often. Screening programs are hpv.routine_screening and hpv.campaign_screening, treatment is hpv.treat_num and hpv.treat_delay, and vaccination is hpv.routine_vx and hpv.campaign_vx. Naming a product as a string, as we do below, picks up the standard version from those files.

Every program reaches women only, unless you pass sex=None for both sexes.

Wire up a screen-and-treat cascade

Real programs are a chain: screen, triage the positives, decide what treatment each woman needs, then deliver it. In HPVsim each link is its own program, and you connect them by naming each one and pointing the next at its outcomes.

import numpy as np
import starsim as ss
import hpvsim as hpv

def screen_and_treat(annual_prob):
    """A VIA screen, triage to a treatment decision, then ablation or excision."""

    screen = hpv.routine_screening(
        product    = 'via',
        prob       = annual_prob,
        start_year = 2015,
        age_range  = [30, 50],
        name       = 'screen',
    )
    assign = hpv.routine_triage(
        product     = 'tx_assigner',
        prob        = 0.9,
        eligibility = lambda sim: sim.interventions['screen'].outcomes['positive'],
        name        = 'assign',
    )
    ablate = hpv.treat_num(
        product     = 'ablation',
        prob        = 0.9,
        eligibility = lambda sim: sim.interventions['assign'].outcomes['ablation'],
        name        = 'ablate',
    )
    excise = hpv.treat_delay(
        product     = 'excision',
        prob        = 0.9,
        delay       = 0.5,                 # Six months from decision to surgery
        eligibility = lambda sim: sim.interventions['assign'].outcomes['excision'],
        name        = 'excise',
    )
    return [screen, assign, ablate, excise]

Three things make this work.

Each link has a name, and the next link’s eligibility is a function that reads the previous link’s outcomes. Positive screens go to assign; the tx_assigner product sorts them into ablation or excision according to how advanced the lesion is; each treatment program picks up its own group.

prob is an annual probability, not a probability per timestep and not a lifetime coverage. prob=0.2 means a woman in the eligible age range has a 20% chance of being screened each year, which averages out to roughly one screen every five years — close to what WHO recommends. prob=0.7 is almost annual screening.

hpv.treat_delay(delay=0.5) puts six months between the treatment decision and the treatment. Women who die or progress in the meantime are never treated, which is what happens when surgical capacity is scarce.

Run it and count the cancers

PARS = dict(
    location       = 'nigeria',
    genotypes      = [16, 18, 'hi5', 'ohr'],
    n_agents       = 2000,
    start          = 1990,
    stop           = 2070,
    ms_agent_ratio = 25,
    rand_seed      = 2,
    verbose        = 0,
)

arms = [
    hpv.Sim(**PARS, label='No screening'),
    hpv.Sim(**PARS, interventions=screen_and_treat(0.2), label='Screen every 5 years'),
    hpv.Sim(**PARS, interventions=screen_and_treat(0.7), label='Screen every year'),
]
msim = ss.MultiSim(arms)
msim.run(verbose=0, parallel=False)
fig = msim.plot('all_hpv_cum_cancers')

Figure(768x576)

baseline = float(msim.sims[0].results.all_hpv.cum_cancers[-1])
for s in msim.sims:
    total = float(s.results.all_hpv.cum_cancers[-1])
    print(f'{s.label:24s} {total:>12,.0f} cancers   '
          f'{(baseline - total) / baseline:>6.1%} averted')

No screening                1,107,529 cancers     0.0% averted
Screen every 5 years          898,669 cancers    18.9% averted
Screen every year             971,483 cancers    12.3% averted

Both screening arms avert cancers, and the more frequent one averts more. Look at how the two averted percentages compare against the difference in screening frequency between them: screening three and a half times as often does not avert three and a half times as much. That diminishing return is the standard finding — most of the benefit comes from screening a woman a few times in her thirties and forties, not from screening her often.

Add therapeutic vaccination

A therapeutic vaccine is offered to women who already have HPV, to help them clear it. HPVsim delivers it either as a campaign to an age group, or to the women a screening program has just found. Here we add the second kind on top of the five-yearly screening program:

def screen_then_txvx(annual_prob):
    cascade = screen_and_treat(annual_prob)
    txvx = hpv.linked_txvx(
        product     = 'txvx1',
        prob        = 0.8,
        eligibility = lambda sim: sim.interventions['screen'].outcomes['positive'],
        name        = 'txvx',
    )
    return cascade + [txvx]

s_screen = hpv.Sim(**PARS, interventions=screen_and_treat(0.2),
                   label='Screen and treat')
s_txvx   = hpv.Sim(**PARS, interventions=screen_then_txvx(0.2),
                   label='Screen, treat, and vaccinate')
msim2 = ss.MultiSim([s_screen, s_txvx])
msim2.run(verbose=0, parallel=False)

for s in msim2.sims:
    print(f'{s.label:30s} {float(s.results.all_hpv.cum_cancers[-1]):>12,.0f} cancers')

Screen and treat                    898,669 cancers
Screen, treat, and vaccinate      1,005,973 cancers

Expect very little difference. The default txvx1 product has an efficacy of 0.01 against precancerous lesions in products_tx.csv — a placeholder standing in for a vaccine that does not yet exist, not a claim about any candidate. If you are modeling a real therapeutic vaccine, put its own efficacy numbers in before reading anything from this comparison.

Build your own product

Every product is a table. To make your own treatment, give one row per health state with the efficacy you want. The states HPVsim knows are susceptible, latent, precin, cin, and cancerous:

import pandas as pd

my_product = pd.DataFrame({
    'name':     'thermal_ablation',
    'state':    ['precin', 'cin', 'cancerous'],
    'genotype': 'all',
    'efficacy': [0.90, 0.90, 0.0],
})
treatment = hpv.tx(df=my_product)

Set genotype to 'all' when efficacy is the same across genotypes, or give one row per genotype when it is not. Diagnostics (hpv.dx) and vaccines (hpv.vx) work the same way; look at the shipped CSVs for the columns each one expects.

The script examples/t05_screen_algorithms.py builds each of the seven screen-and-treat algorithms in the WHO guidelines.