ScenarioPharmacology

Antimalarial drug modelling

There are three different ways of modelling antimalarial drug action in OpenMalaria:

1. A simple success/failure model in which drug treatment either clears all infections in the host (treatSimple in health system / intervention component) or does nothing. This can be extended to include a prophylactic effect which simply removes all new infections for a configurable number of steps in the future.
2. Clearing infections via above method then using Pre-erythrocytic vaccines to probabilistically prevent reinfection (this allows a slightly more realistic prophylactic effect than the above).
3. Explicit pharmacokinetic and pharmacodynamic (PK/PD) modelling of drugs (v33 onwards). A complete xml configuration file using much of this functionality is available.

Explicit pharmacokinetic and pharmacodynamic (PK/PD) models

Required parasite dynamics and health system options

A dynamic model of the course of parasite densities during an infection is required for explicit PK/PD modeling. This is possible only for a limited set of model and health system options as follows:

Within the <model> element the following options must be used (as well as appropriate parameter values):

    <ModelOptions>
      <option name="MOLINEAUX_WITHIN_HOST_MODEL" value="true"/>
      <option name="INCLUDES_PK_PD" value="true"/>
    </ModelOptions>  

Within the health system / intervention component a treatPKPD option must be set e.g. to ensure routine treatment uses a dynamic PKPD model (in this example with drug CQ):

<healthSystem>
          <DecisionTree5Day name="example">
              <pSeekOfficialCareUncomplicated1 value="0.0006"/>
              <pSelfTreatUncomplicated value="0"/>
              <pSeekOfficialCareUncomplicated2 value="0.0006"/>
              <pSeekOfficialCareSevere value="0.48"/>
              <treeUCOfficial>
               <treatPKPD schedule="CQ_3_dose" dosage="CQ"/>
              </treeUCOfficial>
          </DecisionTree5Day>
<healthSystem>

and to carry out mass drug administration (in this example with drug AL) using a dynamic PKPD model the following intervention component is required:

      <component id="MDA">
        <decisionTree>
          <treatPKPD schedule="AL_treat" dosage="AL_dosages"/>
        </decisionTree>
      </component>

For this to function, drug-specific treatment schedules, [dosage tables] (#dosage-tables), PK and PD parameterisations must be specified in the the <pharmacology> element. This has the following structure:

<pharmacology>
  <treatments>
    <schedule name="..."> ... </schedule>
    <schedule name="..."> ... </schedule>
    ...
    <dosages name="..."> ... </dosages>
    <dosages name="..."> ... </dosages>
    ...
  <treatments/>
  <drugs>
    <drug abbrev="..."> ... </drug>
    <drug abbrev="..."> ... </drug>
    ...
  </drugs>
</pharmacology>

Treatment schedules

The treatment schedules determine when drugs are given, relative to the time that treatment is initiated. They are specified by name, and in turn specify the drug, time and "mg" triplets to specify doses. The "mg" value is usually not the total amount of drug given; this is determined in conjunction with a dosage table (the "mg" value may be multiplied, e.g. 1 pill for an infant and 4 for an adult, or multiplied directly by the patient's mass in kg).

	<schedule name="CQ_3_dose"> 
      <medicate drug="CQ" hour="0" mg="10"/>
      <medicate drug="CQ" hour="24" mg="10"/>
      <medicate drug="CQ" hour="48" mg="5"/>
	</schedule>

Drugs are administered according to the specified schedule, but the dosage table may modify the quantity of drug given; for example "CQ_3_dose" may specify three daily administrations of CQ. Poor drug quality and poor adherence may be modelled indirectly by using a dosage table with lower than expected quantities of active ingredient and using a schedule which misses certain doses.

Dosage tables

Dosage tables are also specified by name in the <pharmacology> element. For example, the code below configures one schedule of three doses over three days (using a base amount of active ingredient) and a dosage table multiplying the base quantities in the schedule to get a final quantity (based on body mass in this case, but optionally based on age):

<pharmacology>
  <treatments>
    <schedule name="CQ_3_dose">
      <!-- This is recommend dose per kg (KK) -->
      <medicate drug="CQ" hour="0" mg="10"/>
      <medicate drug="CQ" hour="24" mg="10"/>
      <medicate drug="CQ" hour="48" mg="5"/>
    </schedule>
    <dosages name="CQ">
      <!-- contrived multipliers by DH to give roughly a 1-per-kg multiplier,
      rounding up -->
      <bodymass lowerbound="0" dose_mult="5"/>
      <bodymass lowerbound="5" dose_mult="10"/>
      <bodymass lowerbound="10" dose_mult="20"/>
      <bodymass lowerbound="20" dose_mult="40"/>
      <bodymass lowerbound="40" dose_mult="70"/>
    </dosages>
  </treatments>
  <drugs>
  </drugs>
</pharmacology>

Doses can be determined in one of three ways:

1. According to age. In this case the mg specified in the schedule is a base amount (e.g. one pill or 1/4 pill) and is multiplied according to the relevant age group from the table, e.g. 2: <age lowerbound="5" dose_mult="2"/> (lowerbound is an age in years; multiple age elements should be used).
2. According to weight (mass) in kg. This is the same as for age, except that lowerbound is now a mass in kg and age is replaced with bodymass: <bodymass lowerbound="20" dose_mult="3"/>. Weight (mass) is computed using a function mapping mass onto age in a piece-wise linear manner. (This function does not allow for any variation between hosts or between samples).
3. In terms of mg active ingredient per kg body mass (versions 34 and later). In this case a single multiply element is used: <multiply by="kg"/> and the mg specified in the schedule is mg/kg.

Thus, a dosage tables using age instead of body mass might read:

<dosages name="standard_dosages">
  <age lowerbound="0" dose_mult="5"/>
  <age lowerbound="2" dose_mult="20"/>
  <age lowerbound="6" dose_mult="60"/>
</dosages>

To multiply directly by body-mass (thus making the dose specified in "mg" to be mg/kg), the xml reads:

<dosages name="perkilo">
  <multiply by="kg"/>
</dosages>

Delayed dosing

It is also possible to delay all doses relative to a standard schedule by some number of hours, without creating a new schedule:

<treatPKPD schedule="AL_treat" dosage="AL_dosages" delay_h="84"/>

PK models

The specific PK/PD drug description depends on the structural PK model for the drug. From v34 this can be chosen from any of the following:

1. One-compartment PK model, assuming instantaneous absorption, linear elimination and multiple doses (without lag time) Winter & Hastings (2011).
2. One-compartment absorption-conversion PK model as described in Kay & Hastings (2013). This allows for absorption of a drug from the gut at a user-defined rate, into the unconverted parent form in the serum where it is either eliminated (unconverted) or converted into an active metabolite form at a user-defined rate. The active metabolite is then eliminated at a user-defined rate.
3. Three-compartment PK model with first-order absorption, linear elimination and multiple doses (without lag time) as described by Bertrand & Mentré (2008) (Equation 1.72). The equations were subsequently edited for inclusion into the OpenMalaria code (edited equations ).
4. A two-compartment PK model can be simulated by setting the inter-compartment clearance between the central and peripheral compartment two to zero.

XML and pharmacokinetic/pharmacodynamic (PK/PD) parameter values are available within the pharmacology library as follows:

Drug	One compartment	One Comp. absorption-conversion	Two compartment	Three compartment
artemether	(available)	available and recommended
artesunate	(available)	available and recommended
DHA	available
chloroquine	available
lumefantrine	available
mefloquine	available
piperaquine	(available)		recommended	recommended

In general, the PK parameters of a drug are specified within the element. Each element must contain one element.

<PK>
IN THIS ORDER:
|   <negligible_concentration ... /> 
| EXACTLY ONE OF:
| |   <half_life ... /> 
| | IN THIS ORDER:
| | |   <k ... /> 
| | |   <m_exponent ... /> 
| [ <k_a ... /> ]
| [ <conversion ... /> ]
|   <vol_dist ... /> 
| [ <compartment2 ... /> ]
| [ <compartment3 ... /> ]
</PK>

Some of these parameters are required for only a subset of the four different PK models. More documentation is given below separately for each of the the models. The automatically generated documentation also provides additional details.

One-compartment PK model

A drug with a one-compartment PK model can be specified as follows:

<drug abbrev="MQ">
  <PD>
    <phenotype>
      <max_killing_rate>3.45</max_killing_rate>
      <IC50>0.5</IC50>
      <slope>5</slope>
    </phenotype>
  </PD>
  <PK>
    <negligible_concentration>0.005</negligible_concentration>
    <half_life>13.078</half_life>
    <vol_dist>20.8</vol_dist>
  </PK>
</drug>

Where the PK parameters are:

• negligible_concentration Concentration (mg/l) below which drug's effects are deemed negligible and can be removed from simulation
• half_life time (in days since dosing) at which half the initial concentration remains. Used to calculate elimination rate (which is: ln(2) / half_life). Alternatively, elimination rate can be explicitly specified via k and m_exponent.
• k Constant or distribution associated with elimination rate (k) used to calculate the elimination rate λ, which is calculated as λ = k / (body_mass ^ m_exponent), where body_mass is the patient's weight in kg and m_exponent is the next parameter (an alternative to specifying elimination rate via half_life). The basic form of decay is C(t) = C0 * exp(-λ*t), where t is time in days. If sigma > 0, k is sampled per-human from the log-normal distribution (see below).
• m_exponent Constant used to calculate the elimination rate λ. In the case of a conversion model, this is applied to both the elimination and the conversion rates.
• vol_dist Volume of Distribution in litres/kg. The volume of distribution can be sampled per-human from a log-normal distribution by specifying an additional parameter, sigma (see below).

One-compartment PK model with drug absorption and conversion

Drug absorption and conversion (to an active metabolite) can be incorporated into the one-compartment model by including PK/PD estimates for both the parent drug and active metabolite. To model artemether (AR) with the active metabolite DHA, the PK element would have the following structure:

  <drug abbrev="DHA_AR">
    <PD> ... </PD>
    <PK>
      <k>44.15</k>
      <m_exponent>0.0</m_exponent>
      <vol_dist>15</vol_dist>
    </PK>
  </drug>
  <drug abbrev="AR">
    <PD> ... </PD>
    <PK>
      <k>0.0</k>
      <m_exponent>0.0</m_exponent>
      <k_a>23.98</k_a> <!-- absorption rate -->
      <conversion>
        <metabolite>DHA_AR</metabolite>
        <rate>11.98</rate> <!-- conversion rate -->
        <molRatio>0.9547587</molRatio> <!-- mol. weight ratio -->
      </conversion>
      <vol_dist>46.6</vol_dist>
    </PK>
  </drug>

The parameterisation of the conversion model requires:

• the half-life or k defined as in the simple one compartment model. k may be set to zero for drugs with no direct elimination.

• k_a the absorption rate; this may be large. This also required for two and three compartment models.

• A <conversion> element for the parent drug. This contains:

  • A <metabolite> element which gives the abbreviation of the metabolite (e.g. "DHA" or "DHA_AR"). The metabolite in turn needs to be specified as a separate drug. The specification of the metabolite includes parameters for the half-life or k (elimination of metabolite) (as for the single compartment model), but not k_a (absorption rate). It is not possible for the metabolite to itself undergo conversion with the current models.
  • <rate> is the rate of conversion of parent drug to metabolite.
  • <molRatio> is the ratio of the molecular weight of the metabolite to the molecular weight of the parent.

Two or three-compartment PK models (v34 and later)

To describe a two or three-compartment PK model the <PK>element needs to include estimates of the drugs elimination rate parameter k (see below) and the [m_exponent] (schema-34#elt-m_exponent).

Two new elements, <compartment2> and <compartment3>, are also required to describe the rate (per day) that the drugs move between the PK compartments.

• a12 from central to first peripheral compartment
• a21 from first peripheral compartment to central compartment
• a13 from central to second peripheral compartment
• a31 from second peripheral compartment to central compartment

For example, a two-compartment PK model requires:

<drug abbrev="PPQ"> 
  <PD> ... </PD>
  <PK>
    <negligible_concentration>0.005</negligible_concentration>
    <k mean="0.6242774566473989" sigma="0"/>
    <m_exponent>0.25</m_exponent>
    <vol_dist>173</vol_dist>
    <compartment2>
      <a12 mean="8.46242774566474" sigma="0"/>
      <a21 mean="3.3058064516129035" sigma="0"/>
    </compartment2>
  </PK>
</drug>

while a three-compartment PK model also requires:

<PK>
  <compartment3>
    <a13 mean="43.36156351791531" sigma="0"/>
    <a31 mean="3.8155414012738853" sigma="0"/>
  </compartment3>
</PK>

PD models

The pharmacodynamic (PD) effects of drugs are captured by calculating drug killing factors (1 - parasite survival), which are applied at each time-step to determine how the parasite density at the next step depends on the drug concentration determined from the PK model. In general, immunity impedes parasite growth via an immunity model external to the parasite growth model, and the survival factor calculated from the immunological model is multiplied by the drug survival factor to give an overall impact on parasite survival at each time step.

Combining drug factors from independent drugs

If multiple drugs act independently and simultaneously, each drug is tracked separately. In order to calculate a combined killing factor, the survival factor (one minus the killing factor) is calculated separately for each drug, then these are multiplied to give an overall killing factor.

Combining drug factors in the conversion model

The killing factor is calculated from the parent-metabolite conversion model by taking the maximum killing factor from each drug.

Given factors fP(t) and fM(t) at time t, where a P suffix denotes the parent drug and M the metabolite,

fP(t) = (VP * cP(t) ^ nP) / (cP(t) ^ nP + KnP)
fM(t) = (VM * cM(t) ^ nM) / (cM(t) ^ nM + KnM)

(where V is the max killing rate, c(t) the concentration and Kn = IC50^n).

OpenMalaria calculates the survival factor for some duration as:

exp( -intg[ max(fP(t), fM(t)) ] )

this is slightly different from:

exp( -max( intg[fP(t)], intg[fM(t)] ) )

(the version used by Kay & Hastings (2013) but leads to very similar results and is presumably more computationally efficient because only one numerical integration is required. The duration over which this is calculated may be less than the length of a time step in OpenMalaria

(in the case that drugs are medicated within the time step); in this case the derived factor is taken as the product of the factors for each duration.

Specification of PD parameters

The PD parameters of a drug are specified within the <drug> element. Each <drug> element must contain one element, which in turn specifies one or more <phenotype> elements. In the simple case where spread of drug resistance is not being simulated, only a single <phenotype> element is needed. (When spread of drug resistance is simulated the parasite locus involved in resistance must be indicated and there must be multiple phenotypes identifed by name with different PD parameters. <restriction> elements are then required to link the specific PD parameter set to the associated allele.)

<PD [ locus=string ] >
<phenotype
  [ name=string ]
  >
IN THIS ORDER:
| ( <restriction ... /> )*
|   <max_killing_rate ... /> 
|   <IC50 ... /> 
|   <slope ... /> 
</phenotype>
</PD>

The other PD parameters, which are required, are:

XML element	Symbol	Name	Units	Description
<max_killing_rate>	k1	Maximum killing rate	per day
<IC50>	IC50	IC50	mg/l	Drug concentration at which half the maximum killing rate is achieved. The IC50 can be sampled per-human from a log-normal distribution by specifying an additional parameter, sigma (see below).
<slope>	n	slope		Slope of effect curve

Sampling PKPD parameters from log normal distributions (v34+)

In versions 34 and later of OpenMalaria parameter variability per-human can be added to the following PKPD parameters:

IC50, k, vol_dist, a12, a21, a13, a31.

This is achieved using the sigma (σ) parameter. if the sigma parameter is > 0, the parameter is sampled per-human from the log-normal distribution: ln N( ln(mean) - σ^2 / 2, σ^2). For the parameters k, a12, a21, a13 and a31 the syntax is as follows:

<k mean="0.053" sigma="0.1"/>

where, mean is the untransformed mean and sigma is the standard deviation passed to the log-normal sampler. For the parameters IC50 and vol_dist specification of log-normal sampling requires only the addition of the sigma parameter, e.g.

<IC50="0.2" sigma="0.3"/>

To get an equivalent sample for k in R:

mean <- 0.053
sigma <- 0.1
mu <- log(mean) - sigma*sigma / 2
hist(rlnorm(10000, mu, sigma), breaks=300)

Drug resistance

When explicit PK/PD modeling is implemented, the parasite population can be divided into multiple genotypes with different drug resistance phenotypes, for modeling of the spread of drug resistance and its effects.