Cereal aphid-fungus model

In this section the model is demonstrated step-by-step, each step adding more components (boxes) and complexity (relations between boxes). The plugin section describes the box classes that were developed for this model in particular.

Background

This is a tri-trophic model of the winter wheat-cereal aphid (Sitobion avenae)-fungus (Pandora neoaphidis) system, simulating the dynamics of aphid and fungus from spring until harvest (Saussure et al. 2023). The winter wheat model is rudimentary; the crop is only represented as a growth stage modelled by a regression model fitted to Norwegian growing conditions. Weather data (daily average temperature and maximum relative humidity) drives the model. Optionally, the winter wheat model can be left out, and the crop growth stage instead supplied as an additional column in the weather file.

Below you will find models of increasing complexity, building up to the full model found in the scientific paper. The models are composed of building blocks, implemented as C++ classes. The functionality of model building blocks were implemented in C++. If you really want to get the depth of that, you can download the source code. The classes specific to this model is found in the src/plugins/aphid source code folder.

Each model building block was designed to solve a particular task. The model was constructed from building blocks included in the standard plugin (boxes), together with building blocks developed specifically for this model (in the aphid plugin).

You can run the boxscripts shown in the sub-sections below yourself. For that purpose you must install Universal Simulator. For example, to run the crop_aphid_model.box⏷ script, you write at the Universal Simulator prompt:

> run models/aphid/crop_aphid_model.box

Paste the clipboard at the R prompt to see the output, as described in Running simulations.

Crop model

The crop model can either be weather-driven (by a weather file) or data-driven (by recorded wheat growth stages).

Weather-driven

The weather-driven crop model makes use of the CropGrowthStage and CropIsGrowing models, tailored for Norwegian conditions, together with the generic DayDegrees model. A Calendar and a Records box provides the driving data:

// crop_model_weather_driven.box⏷
Simulation sim {
  Calendar calendar {
    .begin = 01/04/2004
    .end   = 31/08/2004
  }
  Records weather {
    .fileName =  "weather/Aarnes_2004.txt"
  } 
  CropGrowthStage cropGrowthStage {
    .valueAtStart       = 20.0
    .dayDegreesHalfways = 800.0
    .slope              = 2.8
    .dayDegrees = ./time[total]
    CropIsGrowing isGrowing {
      .temperature = weather[Tavg] 
    }
    DayDegrees time{
      .T = weather[Tavg]
      .T0 = 0
      .isTicking = ../isGrowing[value]
    }
  }
  OutputR {
    PageR {
      .xAxis = calendar[date]
      .width = 6
      .height = 3.5
      PlotR {
        .ports = cropGrowthStage[value]
      }
    }
  }
}

Run the model from the Universal Simulator prompt with

> run models/aphid/crop_model_weather_driven.box

and paste to R to see the output:

 

Optionally data-driven

The crop sub-model used in the final, complete model is flexible. It will pick the growth stage from the weather file, if a CropGrowthStage column is available. Otherwise, it uses the value provided by the weather-driven model. This feature was needed for the validation runs of the model, because the field observations were not from Norway. In those cases, the CropGrowthStage model cannot be expected to hold; it is more accurate to use the observed crop growth stages. Here is the boxscript to demonstrate a data-driven crop growth stage:

// crop_model_data_driven.box⏷
Simulation sim {
  Calendar calendar {
    .begin = 01/04/1995
    .end   = 31/08/1995
  }
  Records weather {
    .fileName = "validation/A95-weather-gs.txt"
  } 
  Box wheat{
    CropGrowthStage cropGrowthStageModel {
      .valueAtStart       = 20.0
      .dayDegreesHalfways = 800.0
      .slope              = 2.8
      .dayDegrees = ./time[total]
      CropIsGrowing isGrowing {
        .temperature = weather[Tavg] 
      }
      DayDegrees time{
        .T = weather[Tavg]
        .T0 = 0
        .isTicking = ../isGrowing[value]
      }
    }
    Box cropGrowthStage {
      // Pick growth stage either from the weather file 
         or else from the growth stage model
      &value = if exists(weather[GrowthStage]) then weather[GrowthStage] 
                   else ../cropGrowthStageModel[value]
    }
  }
  OutputR {
    PageR {
      .xAxis = calendar[date]
      .width = 6
      .height = 3.5
      PlotR {
        .ports = cropGrowthStage[value]
      }
    }
  }
}

Run the model from the Universal Simulator prompt with

> run models/aphid/crop_model_data_driven.box

and paste to R to see the output:

Crop-aphid model

In this model, the fungus is left out, so it is basically a model of aphid population dynamics driven bottom-up by the crop. The model is composed of these dedicated model building blocks, which is worth to look up. They each solve one particular task specific to this model:

• CropGrowthStage
• CropIsGrowing
• AphidImmigration
• AphidNetReproduction
• AphidOffspring
• AphidJuvenileSurvival

In addition, you will find these basic building blocks:

• Calendar
• Records
• Stage

After you have loaded the model:

> load models/aphid/crop_aphid_model.box

The list command will show you the model structure:

> list
Simulation sim
  Box param
  Calendar calendar
  Records weather
  Box wheat
    CropGrowthStage cropGrowthStageModel
      CropIsGrowing isGrowing
      DayDegrees time
    Box cropGrowthStage
  Box aphid
    DayDegrees time
    Box density
    AphidImmigration immigration
    AphidNetReproduction netReproduction
    AphidOffspring offspring
    AphidJuvenileSurvival survival
    Box susceptible
      Box nymph
        Stage apterous
        Stage alate
      Box adult
        Stage apterous
          Stage fecundity
        Stage alate
          Stage fecundity
  OutputR 
    PageR 
      PlotR 
    PageR 
      PlotR 
  OutputWriter outputWriter
    OutputSelector selector

The aphid box is just a Box, which means it acts as a simple container of other boxes. Since we are anticipating a model including the fungus soon, we have put the nymph and adult stages inside a susceptible box ('susceptible' means 'uninfected' in epidemiologist speak). Both are further divided into sub-populations of apterous and alate morphs. The adult box holds boxes inside for the individuals (apterous and alate) and their respective reproduction (fecundity).

The param box holds parameter values that must be the same across several boxes, as seen in the boxscript below.

// crop_aphid_model.box⏷
Simulation sim {
  Box param {
    &k = 15
    &juvenileApterousDuration = 172
    &juvenileAlateDuration = 195 
    &adultDuration = 300
    &fecundityDuration = 100
    &fecundity_k = 1
  }
  Calendar calendar {
	...
  }
  Records weather {
    ...
  } 
  Box wheat {
    ...
  }
  Box aphid {
    DayDegrees time{
      .T = weather[Tavg]
      .T0 = 3 
      .Topt = 18 
      .Tmax = 30 
    }  
    Box density {
      &nymphs      = sum(../susceptible/nymph/*[content])
      &adults      = sum(../susceptible/adult/*[content])
      &total       = .[nymphs] + .[adults]
    }
    AphidImmigration immigration {
      .cropGrowthStage = cropGrowthStage[value]
      .toCropGrowthStage = 69
      .immigrationRate = 0.02
      .propExposedImmigrants = 0
      .k = param[k]
    }
    AphidNetReproduction netReproduction {
      .Tmin = 3 
      .Topt = 16.1 
      .Tmax = 30
      .R0opt = 51.6
      .temperature = weather[Tavg]
      .cropGrowthStage = cropGrowthStage[value]
      .optimumCropGrowthStageMin = 59
      .optimumCropGrowthStageMax   = 73
      .optimumCropFactor = 1.6
      .aphidDensity = ../density[total]
      .aphidDensityThreshold = 40 
      .alateFactor = 0.67
    }
    AphidOffspring offspring {
      .offspringTotal = sum(../*/adult/*/fecundity[outflow]) 
      .aphidDensity = ../density[total]
      .cropGrowthStage = cropGrowthStage[value]
    }
    AphidJuvenileSurvival survival{
      .cropGrowthStage = cropGrowthStage[value]
      .temperature = weather[Tavg]
    }
    Box susceptible {
      Box nymph {
        Stage apterous {
          .inflow       = ancestors::*/offspring[apterous]
          .timeStep     = ancestors::*/time[step]
          .growthFactor = ancestors::*/survival[value]
          .k = param[k]
          .duration = param[juvenileApterousDuration]
        }
        Stage alate {
          .inflow = ancestors::*/offspring[alate]
          .timeStep = ancestors::*/time[step]
          .k = param[k]
          .duration = param[juvenileAlateDuration]
          .growthFactor = ancestors::*/survival[value]
        }
      }
      Box adult {
        Stage apterous {
          .inflow = ancestors::*/nymph/apterous[outflow]
          .timeStep = ancestors::*/time[step]
          .duration = param[adultDuration]
          .k = param[k]
          Stage fecundity {
            .inflow = ancestors::*/nymph/apterous[outflow]
            .timeStep = ancestors::*/time[step]
            .duration = param[fecundityDuration]
            .k = param[fecundity_k]
            .growthFactor = ancestors::*/netReproduction[apterous]
          }
        }
        Stage alate {
          .inflow = ancestors::*/immigration[susceptible]
          .timeStep = ancestors::*/time[step]
          .duration = param[adultDuration]
          .k = param[k]
          Stage fecundity {
            .inflow = ancestors::*/immigration[susceptible]
            .timeStep = ancestors::*/time[step]
            .duration = param[fecundityDuration]
            .k = param[fecundity_k]
            .growthFactor = ancestors::*/netReproduction[alate]
          }
        }
      }
    } // susceptible
  } // aphid
  OutputR {
    .width = 5
    .height = 10
    PageR {
      .xAxis = calendar[date]
      PlotR {
        .ports = susceptible/*/*[value] | offspring[output::*]
      }
    }
    PageR {
      .xAxis = cropGrowthStage[value]
      PlotR {
        .ports = susceptible/*/*[value] | offspring[output::*]
        .ggplot = "scale_x_continuous(breaks=10*(0:10))"
      }
    }
  }
}

If you stumble over expressions, such as .[nymph] or ancestors::*/time[step], please read up on port paths.

When you run the model,

> run models/aphid/crop_aphid_model.box

Aphid population dynamics are shown both on a date and a crop growth stage scale. Here are the two plots. First by date:

Then by growth stage:

The alate offspring comes in two waves, the second wave due to decreasing aphid density caused by the ripening of the crop. This may or may not reflect reality. The interaction between aphid density and crop growth stage, which determines aphid fecundity, is complicated.

Crop-aphid model with uncertainty

In this extension of the model above, we include uncertainty on weather and a few model parameters. Every time you run the model, you will get a different output. In this boxscript listing, only the additional boxes and major changes are shown (left-out parts shown as ...):

// crop_aphid_model_ua.box⏷
Simulation sim {
  .iterations = 30
  .silent = TRUE
  SelectFile weatherFiles {
    .folder = "weather"
    .filter = "*.txt"
    .showFileNames = FALSE
  }
  Box random {
    RandomiserStratified randomiser {
    }
    RandomUniformInt k {
      .min = 15
      .max = 30
    }
    RandomUniformInt fileNumber {
      .min = 1
      .max = weatherFiles[numFiles]
    }
    RandomNormal cropAtStart {
      .min = 10
      .max = 30 
    }
    RandomNormal cropHalfways {
      .min = 750 
      .max = 850 
    }
  }
  ...
  Records weather {
    .fileName =  ./files[fileNamePath]   
    .ignoreYear = TRUE
    SelectFile files {
      .folder = "weather"
      .filter = "*.txt"
      .selectFileNumber = random/fileNumber[value]
      .showFileNames = FALSE
    }
  } 
  Box wheat{
    CropGrowthStage cropGrowthStageModel {
      .valueAtStart       = random/cropAtStart[value]
      .dayDegreesHalfways = random/cropHalfways[value]
      .slope              = 2.8
      .dayDegrees = ./time[total]
      CropIsGrowing isGrowing {
        .temperature = weather[Tavg] 
      }
      DayDegrees time{
        .T = weather[Tavg]
        .T0 = 0
        .isTicking = ../isGrowing[value]
      }
    }
    Box cropGrowthStage {
      &value = if exists(weather[GrowthStage]) then weather[GrowthStage] 
                   else ../cropGrowthStageModel[value]
    }
  }
  ...
  OutputR {
    .scripts = "crop_aphid_model_ua.R"
    OutputText {
      .ports = calendar[date] | cropGrowthStage[value] | 
               susceptible/*/*[value] | offspring[output::*] 
    }
  }
}

The random box draws numbers for the uncertain parameters of which there are four:

• k
• fileNumber
• cropAtStart
• cropHalfways

The very first child of random determines the method by which to draw random numbers. In this boxscript, stratified random sampling has been chosen.

The k parameter will achieve a number between 15 and 30 before each iteration of the simulation. It is used by immigration and the four aphid sub-populations. You can use list with the x option to verify this:

> list random/k x

RandomUniformInt k
  >value = 21
    >> sim/aphid/immigration[k]
    >> sim/aphid/susceptible/nymph/apterous[k]
    >> sim/aphid/susceptible/nymph/alate[k]
    >> sim/aphid/susceptible/adult/apterous[k]
    >> sim/aphid/susceptible/adult/alate[k]

The cropGrowthStageModel box uses the random values for cropAtStart and cropHalfways.

fileNumber represents a bit of a hack to choose at random, a file with suffix txt within the weather sub-folder. Two SelectFile boxes do the trick. The first one (weatherFiles) is queried by fileNumber about the total number of candidate files. The second one (files) is used to pick a file by the random number (the files are considered numbered alphabetically), which ultimately is the one opened by weather.

iterations have been set to 30, which means the model will be running 30 times. To avoid an avalanche of status messages on the screen, silent has been set to TRUE. You should try setting it to FALSE (its default value), just to notice the difference.

The OutputR box has been set to run a dedicated R script to process the simulation output:

# crop_aphid_model_ua.R⏷
densities = colnames(sim)[3:9]

M = melt(sim, id.vars=c("iteration", "date"), measure.vars=densities, 
         variable.name="Variable", value.name="Value")
open_plot_window(5,10)
P = ggplot(M, aes(date, Value, colour=Variable, group=iteration)) +
  geom_line(alpha=0.3) +
  labs(y="") +
  theme(legend.position="none") +
  facet_wrap(~Variable, ncol=1, scales="free")
print(P)

M = melt(sim, id.vars=c("iteration", "cropGrowthStage"), measure.vars=densities, 
         variable.name="Variable", value.name="Value")
open_plot_window(5,10)
P = ggplot(M, aes(cropGrowthStage, Value, colour=Variable, group=iteration)) +
  geom_line(alpha=0.3) +
  scale_x_continuous(breaks=10*(0:10)) +
  labs(y="") +
  theme(legend.position="none") +
  facet_wrap(~Variable, ncol=1, scales="free")
print(P)

When you run the simulation,

> run models/aphid/crop_aphid_model_ua.box

the R script will be executed in R, when you paste the clipboard there. The plots are similar to the one before but now show output uncertainty as 30 overlaid curves for each variable:

Aphid dynamics are seen to be less uncertain on a crop growth stage scale than on a date scale.

Crop-aphid-fungus model

We are now ready to tackle the full, tri-trophic model. For that purpose we need to add the fungus. However, it is not represented as a fungus population directly (e.g. as fungal biomass or number of spores in various compartments, such as soil and hosts). Lack of data does not allow us that detail. Instead, the fungus is represented in

• exposed immigrants
• exposed aphids (nymps and adults, apterous and alate)
• aphid cadavers

Remember that exposed means 'infected'. Furthermore, note that we call nymphs, that will develop into alate adults, 'alate'; the correct word is 'alatiform' but we ignore that distinction in the code to make our boxscripts more readable.

Here is the complete boxscript. It has been broken into pieces to allow comments as we go:

// crop_aphid_fungus_model.box⏷
Simulation sim {
  Box param {
    &k = 15
    &juvenileApterousDuration = 172
    &juvenileAlateDuration = 195 
    &adultDuration = 300
    &fecundityDuration = 100
    &fecundity_k = 1
    &lethalTime = 80
  }
  Calendar calendar {
    .begin = 01/04/2004
    .end   = 31/08/2004
  }
  Records weather {
    .fileName = "weather/Aarnes_2004.txt"
  } 
  Box wheat{
    CropGrowthStage cropGrowthStageModel {
      .valueAtStart       = 20.0
      .dayDegreesHalfways = 800.0
      .slope              = 2.8
      .dayDegrees = ./time[total]
      CropIsGrowing isGrowing {
        .temperature = weather[Tavg] 
      }
      DayDegrees time{
        .T = weather[Tavg]
        .T0 = 0
        .isTicking = ../isGrowing[value]
      }
    }
    Box cropGrowthStage {
      &value = if exists(weather[GrowthStage]) then weather[GrowthStage] 
               else ../cropGrowthStageModel[value]
    }
  }
  ...

So far nothing much has changed, only a few additions to param. But here goes! We invoke a Maker box to construct two copies of the box declared in line 5 below. They will be named withoutFungus and withFungus as defined in line 4 and also commented in line 5:

...
Maker system {
  // The system contains two aphid populations, 
  // identical except for the proportion of exposed immigrants
  .names = c("withoutFungus", "withFungus")
  Box { // withoutFungus or withFungus
    DayDegrees time{
      .T = weather[Tavg]
      .T0 = 3 
      .Topt = 18 
      .Tmax = 30 
    } 
    ...

In effect the whole aphid model will be replicated, the two copies having the same structure and sharing all parameter values. They are created as children of system and can referred to by system/withoutFungus and system/withFungus.

We add a fungusTime to keep track of infection development and a few extra ports in aphidDensity to count aphids by their infection status (susceptible, exposed and cadavers). The immigration box asks if the name of its parent is withFungus, in which case 25% of the immigrants will arrive exposed; otherwise zero. Here, name is a BoxScript function. Paths in BoxScript must include a port, but it can be left empty to refer to a box. Hence, we arrive at the expression name(..[]).

...
DayDegrees fungusTime {
  .T = weather[Tavg]
  .T0 = 2 
  .Topt = 18 
  .Tmax = 30 
}
Box aphidDensity {
  &nymphs      = sum(../*/nymph/*[content])
  &adults      = sum(../*/adult/*[content])
  &total       = .[nymphs] + .[adults]
  &susceptible = sum(../susceptible/*/*[content])
  &exposed     = sum(../exposed/*/*[content])
  &cadavers    = sum(../infectious/cadavers[content])
}
AphidImmigration immigration {
  .cropGrowthStage = cropGrowthStage[value]
  .toCropGrowthStage = 69
  .immigrationRate = 0.02
  .propExposedImmigrants = if name(..[]) == "withFungus"
                           then 0.25
                           else 0.0
  .k = param[k]
}
...

The two systems (system/withoutFungus and system/withFungus) do not interact but run in parallel (logically, not in terms of computer processes). This allows comparison between the two systems of the simulation results, since they have everything in common, except for the proportion of exposed immigrants, as defined in lines 19-21 above.

For reproduction and survival, nothing has changed but we need to be carefull with our references. Thus there are two ports on the path aphidDensity[total], namely withoutFungus/aphidDensity[total] and withFungus/aphidDensity[total]. To get the right one, we must use the relative path ../aphidDensity[total]:

...
AphidNetReproduction netReproduction {
  .Tmin = 3 
  .Topt = 16.1 
  .Tmax = 30
  .R0opt = 51.6
  .temperature = weather[Tavg]
  .cropGrowthStage = cropGrowthStage[value]
  .optimumCropGrowthStageMin = 59
  .optimumCropGrowthStageMax   = 73
  .optimumCropFactor = 1.6
  .aphidDensity = ../aphidDensity[total]
  .aphidDensityThreshold = 40 
  .alateFactor = 0.67
}
AphidOffspring offspring {
  .offspringTotal = sum(../*/adult/*/fecundity[outflow]) 
  .aphidDensity = ../aphidDensity[total]
  .cropGrowthStage = cropGrowthStage[value]
}
AphidJuvenileSurvival survival{
  .cropGrowthStage = cropGrowthStage[value]
  .temperature = weather[Tavg]
}
...

The Stage boxes must be extended to allow outflow in a new direction. A certain proportion, set by the model for infectionRate, is leaving susceptible stages to become exposed (further down in the boxscript). First for the nymph stages:

...
Box susceptible {
  Box nymph {
    Stage apterous {
      .inflow       = ancestors::*/offspring[apterous]
      .timeStep     = ancestors::*/time[step]
      .growthFactor = ancestors::*/survival[value]
      .k = param[k]
      .duration = param[juvenileApterousDuration]
      .phaseOutflowProportion = ancestors::*/infectious/infectionRate[value]
    }
    Stage alate {
      .inflow = ancestors::*/offspring[alate]
      .timeStep = ancestors::*/time[step]
      .k = param[k]
      .duration = param[juvenileAlateDuration]
      .growthFactor = ancestors::*/survival[value]
      .phaseOutflowProportion = ancestors::*/infectious/infectionRate[value]
    }
  }
  ...

Then for the adult stages; here, we must remember to let the infected proportion of fecundity follow the adults:

  ...
  Box adult {
    Stage apterous {
      .inflow = ancestors::*/nymph/apterous[outflow]
      .timeStep = ancestors::*/time[step]
      .duration = param[adultDuration]
      .k = param[k]
      .phaseOutflowProportion = ancestors::*/infectious/infectionRate[value]
      Stage fecundity {
        .inflow = ancestors::*/nymph/apterous[outflow]
        .timeStep = ancestors::*/time[step]
        .duration = param[fecundityDuration]
        .k = param[fecundity_k]
        .growthFactor = ancestors::*/netReproduction[apterous]
        .phaseOutflowProportion = ancestors::*/infectious/infectionRate[value]
      }
    }
    Stage alate {
      .inflow = ancestors::*/immigration[susceptible]
      .timeStep = ancestors::*/time[step]
      .duration = param[adultDuration]
      .k = param[k]
      .phaseOutflowProportion = ancestors::*/infectious/infectionRate[value]
      Stage fecundity {
        .inflow = ancestors::*/immigration[susceptible]
        .timeStep = ancestors::*/time[step]
        .duration = param[fecundityDuration]
        .k = param[fecundity_k]
        .growthFactor = ancestors::*/netReproduction[alate]
        .phaseOutflowProportion = ancestors::*/infectious/infectionRate[value]
      }
    }
  }
} // susceptible
...

Notice that the inflow to the nymph stages (apterous and alate) originates from the offspring box, while the inflow to the adult apterous stage originates from individuals finishing their apterous nymph stage. You would expect the same kind of inflow to the adult alate stage; however, alate adults are assumed to fly away immediately. Hence, alate adults originate only has immigrants (which are assumed to stay) received from the immigration box.

No nymphs hatch as exposed (the fungus is only laterally transmitted). Hence the only inflow to the exposed nymph stages is from the infection of susceptible nymphs (just calculated above):

...
Box exposed {
  Box nymph {
    StageAndPhase apterous {
      .timeStep = ancestors::*/time[step]
      .k = param[k]
      .duration = param[juvenileApterousDuration]
      .growthFactor = ancestors::*/survival[value]
      .phaseInflow = ancestors::*/susceptible/nymph/apterous[phaseOutflow]
      .phaseTimeStep = ancestors::*/fungusTime[step]
      .phaseK = param[k]
      .phaseDuration = param[lethalTime]
    }
    StageAndPhase alate {
      .timeStep = ancestors::*/time[step]
      .k = param[k]
      .duration = param[juvenileAlateDuration]
      .growthFactor = ancestors::*/survival[value]
      .phaseInflow = ancestors::*/susceptible/nymph/alate[phaseOutflow]
      .phaseTimeStep = ancestors::*/fungusTime[step]
      .phaseK = param[k]
      .phaseDuration = param[lethalTime]
    }
  } // nymph
  ...

Exposed adults, on the other hand may enter through two pathways: as exposed nymphs developing into exposed adults, or as susceptible adults getting infected. They enter below as the inputs inflow and phaseInflow, respectively:

...
Box adult {
  StageAndPhase apterous {
    .inflow = ancestors::*/nymph/apterous[outflow]
    .timeStep = ancestors::*/time[step]
    .duration = param[adultDuration]
    .k = param[k]
    .phaseInflow = ancestors::*/susceptible/adult/apterous[phaseOutflow]
    .phaseTimeStep = ancestors::*/fungusTime[step]
    .phaseK = param[k]
    .phaseDuration = param[lethalTime]
    StageAndPhase fecundity {
       // No inflow because exposed/nymphs don't reproduce as adults
      .timeStep = ancestors::*/time[step]
      .duration = param[fecundityDuration]
      .k = param[fecundity_k]
      .growthFactor = ancestors::*/netReproduction[apterousExposed] 
      .phaseInflow = ancestors::*/susceptible/adult/apterous/fecundity[phaseOutflow]
      .phaseTimeStep = ancestors::*/fungusTime[step]
      .phaseK = param[k]
      .phaseDuration = param[lethalTime]
    }
  }
  StageAndPhase alate {
    .inflow = ancestors::*/immigration[exposed]
    .timeStep = ancestors::*/time[step]
    .duration = param[adultDuration]
    .k = param[k]
    .phaseInflow = ancestors::*/susceptible/adult/alate[phaseOutflow]
    .phaseTimeStep = ancestors::*/fungusTime[step]
    .phaseK = param[k]
    .phaseDuration = param[lethalTime]
    StageAndPhase fecundity {
      // Exposed immigrants reproduce after arriving
      .inflow = ancestors::*/immigration[exposed]
      .timeStep = ancestors::*/time[step]
      .duration = param[fecundityDuration]
      .k = param[fecundity_k]
      .growthFactor = ancestors::*/netReproduction[alateExposed] 
      .phaseInflow = ancestors::*/susceptible/adult/alate/fecundity[phaseOutflow]
      .phaseTimeStep = ancestors::*/fungusTime[step]
      .phaseK = param[k]
      .phaseDuration = param[lethalTime]
    }
  }
} // adult
...

The exposed box contains one extra child to keep track of cadavers:

  ...
  CadaverConversion succumbed {
    .succumbedApterousNymphs = sum(ancestors::*/nymph/apterous[phaseOutflow])
    .succumbedAlateNymphs    = sum(ancestors::*/nymph/alate[phaseOutflow])
    .succumbedApterousAdults = sum(ancestors::*/adult/apterous[phaseOutflow])
    .succumbedAlateAdults    = sum(ancestors::*/adult/alate[phaseOutflow])
  }
} // exposed
...

Our model in the taxonomy of epidemiological models is an SEI (susceptible-exposed-infectious) model. Hence, to our collection of system boxes (system/susceptible and system/exposed) we now add system/infectious, which contains child boxes to hold cadavers and to calculate infectionRate:

...
Box infectious {
  OnOff isSporulating {
    .x = weather[RHmax]
    .xOn = 95
    .xOff = 999 
  }
  CadaverTime time {
    .isSporulating = ../isSporulating[isOn]
    .timeStep = ancestors::*/fungusTime[step]
    .rhAccelaration = 2
  }
  Stage cadavers {
    .inflow = ancestors::*/exposed/succumbed[cadavers]
    .timeStep = ../time[step]
    .duration = 100
    .k = param[k]
  }
  InfectionRate infectionRate {
    .isSporulating = ../isSporulating[isOn]
    .cadavers = ../cadavers[content]
    .transmissionEfficiency = 0.2
  }
} // infectious
...

The final child of system contains summary statistics, including yield:

      ...
      Box diagnostics {
        Accumulator aphidDays {
          .change = ancestors::*/aphidDensity[total]
        }
        Accumulator cadaverDays {
          .change = ../../infectious/cadavers[content]
        }
        Prevalence prevalence {
          .aphidDensity   = ancestors::*/aphidDensity[total]
          .exposedDensity = ancestors::*/aphidDensity[exposed]
          .cadaverDensity = ancestors::*/aphidDensity[cadavers]
        }
        AphidIndex aphidIndex {
          .nymphs = ancestors::*/aphidDensity[nymphs]
          .adults = ancestors::*/aphidDensity[adults]
        }
        aphid::Yield yield {
          .aphidIndex = ../aphidIndex[value]
          .cropGrowthStage = cropGrowthStage[value]
        }
      } // diagnostics
    } // withoutFungus or withFungus
  } // system
  ...

The two pages of output are both shown as plots in two colums, filled in column-first (that's what direction is for; change it to "row" and notice the difference):

  ...
  OutputR {
    .width = 7
    .height = 7
    PageR {
      .xAxis = calendar[date]
      PlotR {
        .ports = system/*/*/*/Stage::*[value] | Stage::cadavers[value]
        .ncol = 2
        .direction = "col"
        .ggplot = "scale_colour_manual(values=c(rep(red,5), rep(blue,5)))"
      }
    }
    PageR {
      .xAxis = calendar[date]
      PlotR {
        .ports = diagnostics/aphidDays[value] | diagnostics/cadaverDays[value] | 
                 diagnostics/prevalence[output::*] |
                 diagnostics/yield[yield]
        .ncol = 2
        .direction = "col"
        .ggplot = "scale_colour_manual(values=c(rep(red,5), rep(blue,5)))"
      }
    }
  } // OutputR
} // sim (end-of-script)

When we run the model,

run models/aphid/crop_aphid_fungus_model.box

as expected, the fungus is shortening the duration of the aphid outbreak. Cadavers are only present in the system with fungus:

Despite the obvious effect of the fungus on the density of aphids, yield does not seem much improved:

Biocontrol model with uncertainty analysis

The biocontrol_model.box⏷ script adds uncertain parameters, as we already saw in Crop-aphid model with uncertainty, but with even more uncertain parameters, most of them pertaining to the fungus. These are the exact same uncertain 11 parameters that also appear in the scientific paper on the model (). Here is the random box from the biocontrol_model.box⏷ script:

...
Box random {
  RandomiserStratified randomiser {
  }
  RandomUniformInt k {
    .min = 15
    .max = 30
  }
  RandomUniformInt fileNumber {
    .min = 1
    .max = weatherFiles[numFiles]
  }
  RandomNormal cropAtStart {
    .min = 10
    .max = 30 
  }
  RandomNormal cropHalfways {
    .min = 750 
    .max = 850 
  }
  RandomNormal propExposedImmigrants {
    .min = 0.1
    .max = 0.7
  }
  RandomNormal lethalTime {
    .min = 50 
    .max = 115
  }
  RandomNormal immunityCost {
    .min = 0
    .max = 0.4
  }
  RandomNormal sporulationOn {
    .min = 80 
    .max = 99 
  }
  RandomNormal cadaverDuration {
    .min = 48 
    .max = 112 
  }
  RandomNormal timeAcceleration {
    .min = 1
    .max = 3
  }
  RandomNormal transmissionEfficiency {
    .min = 0.05
    .max = 0.5 
  }
}

You can open the boxscript yourself and find how each random parameter is used in the model. The list command with the x option comes in handy (some output abbreviated ...):

> run models/aphid/biocontrol_model.box
...
> list random/RandomBase::* x
RandomUniformInt k
  >value = 30
    >> sim/system/withoutFungus/immigration[k]
    >> sim/system/withoutFungus/susceptible/nymph/apterous[k]
    >> sim/system/withoutFungus/susceptible/nymph/alate[k]
    ...
    >> sim/system/withFungus/immigration[k]
    >> sim/system/withFungus/susceptible/nymph/apterous[k]
    >> sim/system/withFungus/susceptible/nymph/alate[k]
    ...
RandomUniformInt fileNumber
  >value = 6
    >> sim/weather/files[selectFileNumber]
RandomNormal cropAtStart
RandomNormal cropHalfways
RandomNormal propExposedImmigrants
  >value = 0.508888
    >> sim/system/withoutFungus/immigration[propExposedImmigrants]
    >> sim/system/withFungus/immigration[propExposedImmigrants]
RandomNormal lethalTime
  >value = 88.8232
    >> sim/system/withoutFungus/exposed/nymph/apterous[phaseDuration]
    >> sim/system/withoutFungus/exposed/nymph/alate[phaseDuration]
    >> sim/system/withoutFungus/exposed/adult/apterous/fecundity[phaseDuration]
    ...
	>> sim/system/withFungus/exposed/nymph/alate[phaseDuration]
    >> sim/system/withFungus/exposed/adult/apterous/fecundity[phaseDuration]
    >> sim/system/withFungus/exposed/adult/apterous[phaseDuration]
    ...
RandomNormal immunityCost
  >value = 0.0987229
    >> sim/system/withoutFungus/netReproduction[immunityCost]
    >> sim/system/withFungus/netReproduction[immunityCost]
RandomNormal sporulationOn
  >value = 84.5715
    >> sim/system/withoutFungus/infectious/isSporulating[xOn]
    >> sim/system/withFungus/infectious/isSporulating[xOn]
RandomNormal cadaverDuration
  >value = 86.299
    >> sim/system/withoutFungus/infectious/cadavers[duration]
    >> sim/system/withFungus/infectious/cadavers[duration]
RandomNormal timeAcceleration
  >value = 1.90676
    >> sim/system/withoutFungus/infectious/time[rhAccelaration]
    >> sim/system/withFungus/infectious/time[rhAccelaration]
RandomNormal transmissionEfficiency
  >value = 0.305888
    >> sim/system/withoutFungus/infectious/infectionRate[transmissionEfficiency]
    >> sim/system/withFungus/infectious/infectionRate[transmissionEfficiency]
> 

A Biocontrol box has been added near the end of the boxscript, and the output is now sent to a dedicated R script (biocontrol_model.R⏷):

...
Biocontrol biocontrol {
  .aphidPressureWithoutF = withoutFungus/diagnostics/aphidDays[value] 
  .aphidPressureWithF = withFungus/diagnostics/aphidDays[value] 

  .yieldWithoutF = withoutFungus/diagnostics/yield[yield] 
  .yieldWithF = withFungus/diagnostics/yield[yield] 
  
  .cropGrowthStage = cropGrowthStage[value]
  .prevalence = withFungus/diagnostics/prevalence[exposed]
  .cadaverPrevalence = withFungus/diagnostics/prevalence[cadavers]
}
OutputR {
  .scripts = "biocontrol_model.R"
  OutputText {
    .ports = calendar[date] | biocontrol[output::*]
  }
}
...

The biocontrol box compares the outputs of system/withoutFungus and system/withFungus. When you run the model,

> run models/aphid/biocontrol_model.box

it will run for 30 iterations. However, here I have changed it to just one iteration

// biocontrol_model.box⏷
Simulation sim {
  .iterations = 1
  ...

to approach the upcoming complexity in baby steps. Here is the output:

What we see changing with time, as an effect of the fungus, are

• the reduction in aphid pressure (area under the aphid density curve, so-called aphid-days)
• the improvement in yield (in %-points)
• the maximum prevalence of cadavers and the growth stage when that occurred
• the maximum prevalence (of exposed aphids) and the growth stage when that occurred

We also see the percentage of cadavers at growth stage 43, 61 and 73.

Since the model contains uncertain parameters, the course of these nine curves (which we will call model outcomes) will be different for every iteration we run the model. But notice, that the course of the curves is not that important, in fact, we are just interested in the final value for each of the outcomes. They tell us everything (or, at least, they summerise the salient points) about the effects of biocontrol.

When you run the model, you will see 30 curves overlaid for each model outcome:

Imagine now that we picked just the final values for all of these curves. We could then show a histogram for each of the nine outcomes. These histograms would should how uncertain model outcomes are, all due to the uncertainty of the 11 model parameters. If you change iterations to 100, or any number higher then that, you will be rewarded by pastel-coloured histograms:

An obvious question, once you have celebrated your success and posted this beautiful figure (which illustrates an uncertainty analysis) on the wall, is, which of the 11 uncertain parameters are most influential on each of these nine model outcomes? That is the purpose of sensitivity analysis, which we turn to next.

Biocontrol with sensitivity analysis

Setting up the analysis

It will take us only a little editing of the biocontrol_model.box⏷ script to add sensitivity analysis (SA). First, you should realise that SA is very computing intensive, so we need to optimise our sampling strategy for the random numbers. We will shift from stratified random sampling (which optimises the random sampling inside the distribution of each model parameter) to Sobol' quasi-random numbers (which optimise the random sampling inside the 11-dimensional distribution across all 11 model parameters). This is set up at the top of the boxscript, where we replace RandomiserStratified with RandomiserSobolSequence:

// biocontrol_model_sa.box⏷
Simulation sim {
  .iterations = 832 // 832 (for demonstration) or 1703936 (for analysis)
  .silent = TRUE
  .unattended = (sim[iterations] > 1000)
  .stopIteration = !question[accepted]
  .show = (sim[iterations] > 1000)
  PopUp question {
    .title = "You are about to run a simulation that will potentially 
             last many hours (at ~40,000 simulations per hour)"
    .text = "Do you want to continue?"
    .icon = "question"
    .buttons = c("Yes", "No")
  }
  SelectFile weatherFiles {
    .folder = "weather"
    .filter = "*.txt"
    .showFileNames = FALSE
  }
  Box random {
    RandomiserSobolSequence randomiser {
      .doSensitivityAnalysis = TRUE
      .bootstrapSize = 100 // 100 (for demonstration) or 10000 (for analysis)
    }
    RandomUniformInt k {
      .min = 15
      .max = 30
    }
    ...
}

Since your computer is likely to lock up during this lengthy simulation, unattended will be set TRUE if there are more than 1000 iterations. Otherwise, you will get an (unharmful) error message when the simulation has finished, because the clipboard would have been locked as well.

When we are using Sobol' numbers we should be careful to make a balanced sampling of parameter space. To obtain that, the number of simulation iterations must equal \(2^n(p+2)\), where \(n\) is a whole, positive number, and \(p\) is the number of parameters sampled (we've got 11). Hence, these are all valid options:

$$ \begin{split} 2^1\cdot(11+2) &= 26 \\ 2^2\cdot(11+2) &= 52 \\ 2^3\cdot(11+2) &= 104 \\ \ldots \\ 2^{17}\cdot(11+2) &= 1,703,936 \\ \ldots \end{split} $$

The last number listed is the one used in the published paper. That's a weekend job! For a quick demonstration we will stick to 832 iterations but you will also be shown, what came out of the weekend job. If you set iterations to an invalid number, the simulation will stop immediately if you try to run it. As a service, you will get suggestions to valid values for iterations near the number you picked. So, if you want in the range of, say, 50,000 iterations, write that, and you will be told nearby valid numbers.

To carry out the sensitivity analysis on the simulation output, we must set doSensitivityAnalysis to TRUE and choose a reasonable bootstrapSize to do statistics on the sensitivity indices. A bootstrapSize of 10,000 seems to be a standard choice in literature but here we'll use 100 for demonstration.

The Popup box is there to warn the user and give the option of regretting ever having started the simulation. The answer is used by sim to stop the simulation in the first step, if stopIteration is TRUE.

Showing the results

The OutputR box contains four child boxes. Here are the first two:

...
OutputR {
  .saveDataFrame = TRUE
  OutputSelector {
    .final = TRUE
  }
  PageR {
    .xAxis = random/*[value]
    PlotR {
      .ports = biocontrol[output::*]
      .maxData = 1000
      .ggplot = "geom_smooth(colour='yellow')"
    }
  }
  ...
}

We ask for the output to be saved in a data frame (saveDataFrame is TRUE) to keep it for later; the simulation may have taken many hours and we don't want to lose the output. We invoke the OutputSelector to write only the final values of each iteration. We are going to use the outputs generated by the Biocontrol box, and for those the final values nicely summarise the effect of biocontrol.

The first PageR box puts all 11 uncertain parameters on the x-axis and all 9 biocontrol outputs on the y-axis. Let's decipher those two paths for xAxis and ports. But first we need to run the model (this took 20 seconds on my computer),

> run models/aphid/biocontrol_model_sa.box

We can use the find command to find out what's on the path specified for the xAxis above (random/*[value]):

> find random/*[value]
Port int    sim/random/k[value]
Port int    sim/random/fileNumber[value]
Port double sim/random/cropAtStart[value]
Port double sim/random/cropHalfways[value]
Port double sim/random/propExposedImmigrants[value]
Port double sim/random/lethalTime[value]
Port double sim/random/immunityCost[value]
Port double sim/random/sporulationOn[value]
Port double sim/random/cadaverDuration[value]
Port double sim/random/timeAcceleration[value]
Port double sim/random/transmissionEfficiency[value]

Check. Those are the 11 parameters we want on the x-axis. Now for biocontrol[output::*]:

> find biocontrol[output::*]
Port double sim/biocontrol[aphidPressureDifference]
Port double sim/biocontrol[yieldImprovement]
Port double sim/biocontrol[percentageCadaversGs43]
Port double sim/biocontrol[percentageCadaversGs61]
Port double sim/biocontrol[percentageCadaversGs73]
Port double sim/biocontrol[maxCadaverPrevalence]
Port double sim/biocontrol[maxCadaverPrevalenceGS]
Port double sim/biocontrol[maxPrevalence]
Port double sim/biocontrol[maxPrevalenceGS]

Check. We've got nine of those.

If you were running the simulation unattended (see above), you must type in the clip command to fill the clipboard after the simulation has finished. When you paste that into R, R will do the bootstrapping on the sensitivity indices. This took 15 seconds on my machine. The first thing you will see, is this impressive 11 by 9 plot:

Notice that in the PlotR box, we set maxData to 1000 to show no more than the first 1000 lines of output; we don't want thousands of dots in those tiny plots. The R code given to the ggplot input is added to the plot in R (where it is produced by R's ggplot function), which is why the yellow trend lines appear on top of the dots. For a polished version of this plot, see the published Figure 6.

The other two plots are defined like this:

PageR {
  .xAxis = random/*[value]
  PlotR {
    .ports = biocontrol[output::*]
    .type = "SobolConvergence"
  }
}
PageR {
  .xAxis = random/*[value]
  PlotR {
    .ports = biocontrol[output::*]
    .type = "SobolIndices"
  }
}

What sets them apart is their type. These are specialised plots meant to be used only in connection with a sensitivity analysis. For these plots, you will need plenty of iterations to give a correct picture of model sensitivities. Here, we are showing how the plots look after 1,703,936 iterations.

The first plot is a convergence plot:

You use it to judge whether the two Sobol' sensitivity indices (total and first-order) have stabilised, or if you need to increase the number of iterations. The number on the x-axis is the so-called sample size \((N)\) of the analysis. It is \(N=2^n\) defined by the number of iterations = \(2^n(p+2)\) (see Setting up the analysis). The first column shows the sums across all 11 parameters for both indices. They are still changing ever so slightly going from \(N=2^{15}\) to \(N=2^{16}\) for some of the nine outcome variables (arranged in rows). We judge that \(N=2^{17}=131,072\) is a sufficient sample size for our analysis. You can find a more focused plot in the published Figure 8.

The second plot shows the total and first-order sensitivity indices for each of the nine outcome variables:

The 11 parameters have been sorted for each of the nine outcome variable to show the most influential ones at the top. The error bars show the 95% confidence limits estimated by bootstrapping. Again, a more focused plot was published as Figure 5. You can tell a long story based upon this sensitivity analysis. For this we refer to our paper (Saussure et al. 2023).

Reproduce published figures

Instructions

Here you will find instructions how to reproduce all the figures published in Saussure et al. (2023).

Some of the figures rely on the sim or S data frames, generated by the big (1,703,936 iterations) sensitivity analysis. You can generate those two data frames yourself, using the biocontrol_model_sa.box script with iterations changed to 1703936 (see Setting up the analysis). Just let your computer work it out over a weekend. Alternatively, you can download them (download sim | download S). Either way, you will end up with two files with suffix (file type) .Rdata . Put those two files in a folder of your choice.

If an .Rdata file is needed to produce a figure, it will be noted below, together with the name of the R script producing the figure. In that R script, you must change the file path at the top to make it read your .Rdata file. Subsequently run the script in R to generate the figure. As an example, below are the steps to generate Figure 3.

First, find the location of the input folder using the get folders command:

> get folders

and open the figure_3.R⏷ script found there (or download it through the provided ⏷ link).

You will have to change three variables, all clearly marked at the top of the R script, designating

• the file name and path leading to your sim or S data frame, possibly both
• the name and path leading to the input/scripts/begin.R script, included with your installation of Universal Simulator
• the path to the folder where you want the graphics files to be written (the figures will be shown both on the screen in R and written to files)

Here is an example:

# Load your sim data
sim_data_file = "~/ipm/aphid-biocontrol-sim.Rdata"

# Load standard script
source("C:/users/joe_cool/UniversalSimulatorHome/input/scripts/begin.R")

# Output folder
output_folder = "~/ipm/figures"

Note: The tilde ~ is a shorthand for your Documents folder. It works even on Windows!

Now save and run the R script.

Figure 2

You create this figure by running a script that is nearly identical to the biocontrol_model.box⏷ script. Here goes

> run models/aphid/figure_2.box

The script creates four figures, two on the screen (one in colours, one in grey tones) and two in graphics files (one PNG and one EPS). The figure_2.R⏷ script creates the figures, writing the graphics files to your Universal Simulator output folder; the precise location is reported by the R script, so that you can easily retrieve them. Here is the PNG file:

Since the boxscript contains uncertain parameter values, you will get different results every time you run the script.

Figure 3

This figure is generated by the figure_3.R⏷ script. Follow the instructions in Reproducing figures from saved data and watch the result:

Figure 4

This figure is generated by the figure_4.R⏷ script. Follow the instructions in Reproducing figures from saved data and watch the result:

Figure 5

This figure is generated by the figure_5.R⏷ script. Follow the instructions in Reproducing figures from saved data and watch the result:

Hey, wait. That didn't quite work out like that. The figure you generated is missing all the labels. That's because we handicrafted them afterwards, as we found it impossible to get these nice-looking labels in R. For proper labelling, you will find a text file with all figure data in fig-5-bw-table.txt, which was written to the output_folder by the figure_5.R⏷ script.

The figure depicts Sobol’ indices \((S_{Ti}\): dark grey; \(S_i\) : light grey\()\) showing model sensitivity to the nine uncertain model parameters, in terms of three model outcomes: (a) peak prevalence of exposed aphids, (b) peak prevalence of cadavers, and (c) yield improvement due to biocontrol. The model parameters are listed in order of importance for each outcome according to \(S_{Ti}\).

Figure 6

This figure is generated by the figure-6.R⏷ script. You need to have the mgcv package installed in R before you run the script. Otherwise, follow the instructions in Reproducing figures from saved data and watch the result:

If you study the figure_6.R⏷ script, you will find that many a trick was pulled to produce this figure! The mathematical symbols were added afterwards outside R.

The figure depicts the response of three model outcomes to variation in the seven most influential model parameters (excluding weather).

Figure 7

This figure is generated by the figure_7.R⏷ script. Follow the instructions in Reproducing figures from saved data and watch the result:

Figure 8

This figure is generated by the figure_8.R⏷ script. Follow the instructions in Reproducing figures from saved data and watch the result:

Figure 9

This figure is generated by the figure_9.R⏷ script. Follow the instructions in Reproducing figures from saved data and watch the result:

Figure 10

You create this figure by running the boxscript,

> run models/aphid/figure_10.box

to get the result:

Figure 11

You create this figure by running the boxscript,

> run models/aphid/figure_11.box

to get the result:

Aphid plugin

The aphid plugin contains the model building blocks needed to run the cereal aphid-fungus model. The functionality of the building blocks (classes) is described below

You can run the boxscripts shown in the sub-sections below yourself. For that purpose you must install Universal Simulator. For example, to run the aphid_juvenile_survival.box⏷ script, you write at the Universal Simulator prompt:

> run models/aphid/aphid_juvenile_survival.box

Paste the clipboard at the R prompt to see the output as described in Running simulations.

AphidImmigration

Interface

Inputs	Type	Default	Purpose / Expression
cropGrowthStage	double	0.0 Zadoks	Current crop growth stage
fromCropGrowthStage	double	31.0 Zadoks	Immigration begins as this growth stage
toCropGrowthStage	double	80.0 Zadoks	Immigration ends as this growth stage
immigrationRate	double	0.0 per tiller d-1	Immigration rate in the season
propExposedImmigrants	double	0.0	Proportion of infected immigrants [0;1]
k	int	0 positive int	Has to match k of the receiving StageAndPhase box
Outputs
total	double	0.0 per tiller d-1	Total immigration rate
susceptible	double	0.0 per tiller d-1	Immigration rate of susceptible aphids
exposed	vec_double	c() per tiller d-1	Immigration rate of exposed aphids

Usage

The AphidImmigration box computes the daily immigration rate of susceptible and exposed aphids (i.e., aphids uninfected or infected by the fungus).

Immigration will be ongoing when cropGrowthStage is between fromCropGrowthStage and toCropGrowthStage at a rate of immigrationRate. The proportion propExposedImmigrants will be considered exposed.

Since the exposed aphids have an age structure reflecting the age of the infection, the exposed output from AphidImmigration needs to be a vector, while the susceptible output is a simple number. The length of the exposed vector is set by the k input.

This boxscript contains a CropGrowthStage model to drive the AphidImmigration model. The default fromCropGrowthStage=31 is maintained while fromCropGrowthStage is set to 69:

// aphid_immigration.box⏷
Simulation sim {

  Calendar calendar {
    .begin = 01/04/2004
    .end   = 31/08/2004
  }

  Records weather {
    .fileName =  "weather/Aarnes_2004.txt"
  } 

  Box wheat{
    CropGrowthStage cropGrowthStage {
      .valueAtStart       = 20.0
      .dayDegreesHalfways = 800.0
      .slope              = 2.8
      .dayDegrees = ./time[total]
      CropIsGrowing isGrowing {
        .temperature = weather[Tavg] 
      }
      DayDegrees time{
        .T = weather[Tavg]
        .T0 = 0
        .isTicking = ../isGrowing[value]
      }
    }
  }
  
  Box aphid {
    AphidImmigration immigration {
      .cropGrowthStage = cropGrowthStage[value]
      .toCropGrowthStage = 69
      .immigrationRate = 0.02
      .k = 15
    }
  }
  
  OutputR {
    PageR {
      .xAxis = calendar[date]
      &exposed = sum(aphid/immigration[exposed])
      PlotR {
        .ports = wheat/cropGrowthStage[value] | 
                 aphid/immigration[susceptible] | ..[exposed]
      }
    }
  }
}

Run the boxscript with

> run models/aphid/aphid-immigration.box

and the output will show the expected pattern, keeping the immigration of exposed aphids at zero, since propExposedImmigrants kept its default value of zero:

AphidIndex

Interface

Inputs	Type	Default	Purpose / Expression
nymphs	double	0.0 per tiller	Aphid nymph density \((N_{nymph})\)
adults	double	0.0 per tiller	Aphid adult density \((N_{adult})\)
Outputs
value	double	0.0	Index value \(y\)

Background

The AphidIndex expresses the severity of the aphid attack by the index defined by Wratten et al. 1979,

$$ y=log_{10}(\,0.5*N_{nymph}+N_{adult}+0.01) $$

It is the logarithm of aphid density with nymphs counting half of the adults. The index is -2 at zero aphid density. It is used as an input to calculate Yield.

Usage

The usage of AphidIndex is demonstrated in the crop-aphid-fungus model.

AphidJuvenileSurvival

Interface

Inputs	Type	Default	Purpose / Expression
cropGrowthStage	double	0.0 Zadoks	Crop growth stage \((G)\)
temperature	double	0.0 oC	Daily average temperature \((T)\)
Outputs
value	double	0.0	Juvenile survival \((y)\) [0;1] d-1

Background

The AphidJuvenileSurvival model computes the daily survival rate of nymphs from the inputs cropGrowthStage and temperature:

$$ y=\left\{ \begin{matrix} 0.944-3.32\cdot10^{-10}\,exp(0.726\cdot T) &\text{for}&G<73\\ 0.45&\text{for}&73\le G<80 \\ 0&\text{for}&80\le G \end{matrix} \right. $$

The model parameters were estimated by Duffy et al. (2017).

Usage

You use AphidJuvenileSurvival as a building block in an aphid population dynamics model. This boxscript runs through temperature from 15 to 35℃ and computes the survival at growth stage 60:

// aphid_juvenile_survival.box⏷
Simulation sim {
  .steps = temperature[steps]
  SequenceByStep temperature {
    .min = 15
    .max = 35
    .by = 0.5
  }
  AphidJuvenileSurvival survival {
    .cropGrowthStage = 60
    .temperature = temperature[value]
  }
  OutputR {
    PageR {
      .xAxis = temperature[value]
      PlotR {
        .ports = survival[value]
      }
    }
  }
}

Run the boxscript with

> run models/aphid/aphid-juvenile-survival.box

and you will get the output:

AphidNetReproduction

Interface

Inputs	Type	Default	Purpose / Expression
R0opt	double	51.6 per capita	Optimum net reproduction \((R_0^{opt})\)
Tmin	double	3.0 oC	Minimum temperature under which no reproduction occur \((T_{min})\)
Tmax	double	30.0 oC	Maximum temperature over which no reproduction occur anymore \((T_{max})\)
Topt	double	16.1 oC	Optimum temperature for reproduction \((T_{opt})\)
temperature	double	0.0 oC	Daily average temperature \((T)\)
cropGrowthStage	double	0.0 Zadoks	Crop growth stage \((G)\)
optimumCropGrowthStageMin	double	0.0 Zadoks	The crop is optimal for reproduction from this growth stage \((G_{min})\)
optimumCropGrowthStageMax	double	0.0 Zadoks	The crop is optimal for reproduction until this growth stage \((G_{max})\)
optimumCropFactor	double	0.0 unitless	Fecundity increased by this factor when crop is optimal \((c_{crop})\)
alateFactor	double	0.0 unitless	Factor to correct alate relative to apterous fecundity \((c_{alate})\)
aphidDensity	double	0.0 per tiller	Aphid density \((N)\)
aphidDensityThreshold	double	0.0 per tiller	Density threshold when net reproduction is zero \((N_{max})\)
immunityCost	double	0.0	Relative reduction in reproduction when exposed \((\nu)\) [0;1]
Outputs
apterous	double	0.0 per capita	Net reproduction for susceptible apterous aphids \((R_0^{aptS})\)
alate	double	0.0 per capita	Net reproduction for susceptible alate aphids \((R_0^{alaS})\)
apterousExposed	double	0.0 per capita	Net reproduction for infected apterous aphids \((R_0^{aptE})\)
alateExposed	double	0.0 per capita	Net reproduction for infected alate aphids \((R_0^{alaE})\)

Background

Aphid fecundity depends on temperature, crop growth stage and aphid density. The AphidNetReproduction model calculates life time fecundity \(R_0\) for a female according to morph (apterous or alate) and infection status (susceptible or exposed).

The temperature response is modelled as a bi-linear function,

$$ f(T) = \left\{ \begin{matrix} \begin{matrix} \frac{R_0^{opt}}{T_{opt}-T_{min}}T - \frac{R_0^{opt}\,T_{min}}{T_{opt}-T_{min}} &\text{for}&T_{min}\le T\lt T_{opt}\\ \frac{R_0^{opt}}{T_{opt}-T_{max}}T - \frac{R_0^{opt}\,T_{max}}{T_{opt}-T_{max}}&\text{for}&T_{opt}\le T\lt T_{max} \\ 0&\text{for}&T\lt T_{min} \vee T_{max}\le T \end{matrix} \end{matrix} \right\} $$

Density-dependence is simply linear up until the threshold,

$$ g(N)=\left\{ \begin{matrix} 1-N/N_{max}&\text{for}&N \lt N_{max} \\ 0&\text{for}&N\ge N_{max} \end{matrix} \right. $$

The effect of crop growth stage is to increase fecundity with the optimal growth stage range:

$$ h(G)=\left\{ \begin{matrix} 1&\text{for}&G\lt G_{min}\\ c_{crop}&\text{for}&G_{min}\le G\lt G_{max} \\ 1&\text{for}&G_{max}\le G<80 \\ 0&\text{for}&G\ge 80 \end{matrix} \right. $$

To achieve \(R_0\) for the four combinations of mother morph and infection status, the effects of temperature \(f(T)\), density \(g(N)\) and crop \(h(G)\) were multiplied together with the effect of morph \(c_{alate}\) and immunity cost \(\nu\):

$$ \begin{split} R_0^{aptS} &= f(T)\,g(N)\,h(G) \\ R_0^{alaS} &= f(T)\,g(N)\,h(G)\,c_{alate} \\ R_0^{aptE} &= f(T)\,g(N)\,h(G)\,(1-\nu) \\ R_0^{alaE} &= f(T)\,g(N)\,h(G)\,c_{alate}\,(1-\nu) \end{split} $$

This boxscript demonstrates the functionality of AphidNetReproduction at different temperatures:

// aphid_net_reproduction.box⏷
Simulation sim {
  .steps = temperature[steps]
  SequenceByStep temperature {
    .min = 0
    .max = 35
    .by = 1
  } 
  AphidNetReproduction reproduction {
    .temperature = temperature[value]
    .Tmin = 3 
    .Topt = 15
    .Tmax = 30
    .R0opt = 50
    .cropGrowthStage = 60
    .optimumCropGrowthStageFrom = 59
    .optimumCropGrowthStageTo   = 73
    .optimumCropFactor = 1.6
    .aphidDensity = 10
    .aphidDensityThreshold = 40 
    .alateFactor = 0.67
    .immunityCost = 0.2
  }
  OutputR {
    PageR {
      .xAxis = temperature[value]
      PlotR {
        .ports = reproduction[output::*]
        .layout = "merged"
      }
    }
  }
}

Run the boxscript with

> run models/aphid/net-reproduction.box

and you will get the output:

We can check the calculation at optimum temperature with the parameter values set in the boxscript:

$$ \begin{split} R_0^{apt} &= 50\cdot(1-10/40)\cdot1.6=60\\ R_0^{ala} &= 50\cdot(1-10/40)\cdot1.6\cdot0.67=40.2\\ R_0^{aptE} &= 50\cdot(1-10/40)\cdot1.6\cdot(1-0.2)=48\\ R_0^{alaE} &= 50\cdot(1-10/40)\cdot1.6\cdot0.6\cdot(1-0.2)=32.16 \end{split} $$

Usage

The usage of AphidNetReproduction is demonstrated in the crop-aphid model.

AphidOffspring

Interface

Inputs	Type	Default	Purpose / Expression
offspringTotal	double	0.0 per tiller	Total no. of offspring produced by susceptible aphids \((n_{total})\)
aphidDensity	double	0.0 per tiller	Aphid density \((N)\)
cropGrowthStage	double	0.0 Zadoks	Crop growth stage \((G)\)
Outputs
apterous	double	0.0 per tiller	Total no. of apterous offspring produced \((n_{apt})\)
alate	double	0.0 per tiller	Total no. of alate offspring produced \((n_{ala})\)
alateProportion	double	0.0	Proportion of alate offspring \((p_{ala})\) [0;1]

Background

The AphidOffspring model splits offspringTotal into apterous and alate (strictly, 'alatiform') offspring in response to aphidDensity and cropGrowthStage based on the empirical relation Watt and Dixon, 1981:

$$ \begin{align*} p_{ala} &= \left[ (2.603N + 0.847G - 27.189)/100 \right]_0^1\\ n_{apt} &= (1-p_{ala})\,n_{total} \\ n_{ala} &= p_{ala}\,n_{total} \end{align*} $$

where \([\ldots]_0^1\) denotes that the expression is limited to the closed interval \([0;1]\).

Usage

The usage of AphidOffspring is demonstrated in the crop-aphid model.

Biocontrol

Interface

Inputs	Type	Default	Purpose / Expression
aphidPressureWithoutF	double	0.0 aphid days	Accumulated aphid pressure without fungus
aphidPressureWithF	double	0.0 aphid days	Accumulated aphid pressure with fungus
yieldWithoutF	double	0.0	Relative yield without fungus [0;1]
yieldWithF	double	0.0	Relative yield witt fungus [0;1]
cropGrowthStage	double	0.0 Zadoks	Current crop growth stage
prevalence	double	0.0 %	Prevalence of exposed aphids
cadaverPrevalence	double	0.0 %	Prevalence of cadavers
Outputs
aphidPressureDifference	double	0.0 aphid days	Difference in aphid pressure caused by fungus
yieldImprovement	double	0.0 %-points	Improvement in yield when controlled
percentageCadaversGs43	double	0.0 %	Percentage cadavers at GS 43
percentageCadaversGs61	double	0.0 %	Percentage cadavers at GS 61
percentageCadaversGs73	double	0.0 %	Percentage cadavers at GS 73
maxCadaverPrevalence	double	0.0 %	Peak cadaver prevalence before GS80
maxCadaverPrevalenceGS	double	0.0 Zadoks	Crop growth stage at peak cadaver prevalence before GS80
maxPrevalence	double	0.0 %	Peak prevalence before GS80
maxPrevalenceGS	double	0.0 Zadoks	Crop growth stage at peak prevalence before GS80

Background

The Biocontrol box summarises the effects of the fungus in terms of aphid pressure, yield and the prevalence of cadavers and exposed aphids.

Usage

The usage of Biocontrol is demonstrated in the biocontrol model.

CadaverConversion

Interface

Inputs	Type	Default	Purpose / Expression
succumbedApterousNymphs	double	0.0 per tiller	New apterous nymph cadavers
succumbedAlateNymphs	double	0.0 per tiller	New alate nymph cadavers
succumbedApterousAdults	double	0.0 per tiller	New apterous adult cadavers
succumbedAlateAdults	double	0.0 per tiller	New alate adult cadavers
Outputs
cadavers	double	0.0 per tiller	Cadavers in standardised units
count	double	0.0 per tiller	Number of cadavers

Background

The inputs are newly (i.e, during the latest time step) produced cadavers of different life stages and morphs. The outputs are the total number of cadavers \(C_n\) and the total number, standardised to adult cadavers \(C_{std}\):

$$ \begin{align*} C_{n} &= N_{apt} + N_{ala} + A_{ala} + A_{apt} \\ C_{std} &= 0.5\left( N_{apt} + N_{ala} \right) + 0.8 A_{ala} + A_{apt} \end{align*} $$

The rationale behind \(C_{std}\) is that alates produce less (20%; cf. Hemmati et al., 2001) spores than apterae and that nymphs produce fewer (50%; guessed) spores than adults.

Usage

The usage of CadaverConversion is demonstrated in the crop-aphid-fungus model.

CadaverTime

Interface

Inputs	Type	Default	Purpose / Expression
isSporulating	bool	FALSE boolean	Are cadavers sporulating?
timeStep	double	0.0 DD	Time step \((\tau)\)
rhAccelaration	double	0.0 -	Acceleration factor above RH threshold \((h)\)
Outputs
step	double	0.0 DD	RH-corrected time step \((\tau_{corrected})\)
total	double	0.0 DD	Accumulated RH-corrected time steps since reset

Background

Cadavers degrade on a day-degree scale. When they are sporulating this degradation is accelerated by a factor rhAccelaration. The CadaverTime box simple multiplies timeStep with rhAccelaration if isSporulating is true:

$$ {\tau}_{corrected} = \left\{ \begin{matrix} h\tau \,\,&\text{if is sporulating} \\ \tau \,\,&\text{otherwise} \end{matrix} \right. $$

Usage

The usage of CadaverTime is demonstrated in the crop-aphid-fungus model.

CropGrowthStage

Interface

Inputs	Type	Default	Purpose / Expression
valueAtStart	double	20.0 Zadoks	Growth stage at the beginning of the growth season \((G_0)\)
dayDegrees	double	0.0 DD	Day-degrees passed since growth season started \((\tau)\)
dayDegreesHalfways	double	720.0 DD	Time when development is half completed \((\tau_{50})\)
slope	double	2.8 DD-1	Max. development rate \((g)\)
Outputs
value	double	0.0 Zadoks	Crop growth stage \((G)\)

Background

The CropGrowthStage is an empirical model driven by day-degrees and fitted to Norwegian data. It follows the sigmoid curve of the log-logistic function

$$ G = G_0 + \frac{G_{max}-G_0}{1+exp\{1+g[\,ln(\tau)-ln(\tau_{50})]\,\}} $$

with \(G_{max}=99\).

Usage

The usage of CropGrowthStage is demonstrated together with AphidImmigration.

CropIsGrowing

Interface

Inputs	Type	Default	Purpose / Expression
temperature	double	0.0 oC	Daily average temperature
T0	double	5.0 oC	Threshold that triggers crop growth
Outputs
value	bool	FALSE	Is the crop currently growing?

Background

The crop goes out of hibernation to continue growth and development, when experiencing five consecutive days with temperature above the threshold T0. The output value is false at the beginning of the simulation and remains true onwards, once the condition has been fulfilled.

Usage

The usage of CropIsGrowing is demonstrated in the crop model.

InfectionRate

Interface

Inputs	Type	Default	Purpose / Expression
isSporulating	bool	FALSE boolean	Are cadavers sporulating?
cadavers	double	0.0 per tiller	Sporulating cadavers \((C)\)
transmissionEfficiency	double	0.0 per cadaver per day	Transmission efficiency \((\epsilon)\)
Outputs
value	double	0.0	Proportion of hosts getting exposed in one day \((\epsilon_{finite})\) [0;1]

Background

Infection spreads at the instantaneous rate $\epsilon$ per day. This is transformed into a finite rate per day applying the functional response of Nicholson and Bailey 1935:

$$ {\epsilon}_{finite} = \left\{ \begin{matrix} 1-exp\left(-\epsilon C\right) \,\,&\text{if is sporulating} \\ 0 \,\,&\text{otherwise} \end{matrix} \right. $$

Usage

The usage of InfectionRate is demonstrated in the crop-aphid-fungus model.

Prevalence

Interface

Inputs	Type	Default	Purpose / Expression
aphidDensity	double	0.0 per tiller	Current total density of live aphids
exposedDensity	double	0.0 per tiller	Current density of exposed aphids
cadaverDensity	double	0.0 per tiller	Current density of cadavers
Outputs
exposed	double	0.0 %	Percentage exposed aphids out of all live aphids
cadavers	double	0.0 %	Percentage cadavers out of aphids+cadavers

Background

The Prevalence box calculates the current percentages of exposed aphids and of cadavers.

Usage

The usage of Prevalence is demonstrated in the crop-aphid-fungus model.

Yield

Interface

Inputs	Type	Default	Purpose / Expression
cropGrowthStage	double	0.0 Zadoks	Crop growth stage \((G)\)
aphidIndex	double	0.0	Aphid density index of Entwistle-Dixon-Wratten \((I)\)
Outputs
yield	double	0.0	Yield relative to uninfested crop \((Y)\) [0;1]
loss	double	0.0	Proportional yield loss due to aphids \((Y_L)\) [0;1]
lossRate	double	0.0	Daily loss rate \((y)\) [0;1]

Background

The Yield class implements the yield loss model of Entwistle and Dixon 1987, which takes cropGrowthStage (from CropGrowthStage) and aphidIndex (from AphidIndex) as inputs. The crop stage \(G\) is used in a lookup table to find the coefficient \(E\):

\(G\)	\(E\)
[0; 53)	0.075
[53; 69)	0.205
[69; 71)	0.075
[71; 73)	0.056
[73; 77)	0.037
[77;99]	0.012

The daily loss rate \(r \in [0;1]\) is found by multiplying the aphid index \(A\) by \(E\):

$$ r= \left\{ \begin{matrix} AE/100 &\text{for}&A>0 \\ 0 &\text{for}&A\le 0 \end{matrix} \right. $$

The Yield object keeps track of the accumulated loss by updating the relative yield \(Y \in [0;1]\) and loss \(Y_L \in [0;1]\) in every time step (1 day):

$$ \begin{align*} Y &\gets (1-r)Y \\ Y_L &= 1 - Y \end{align*} $$

Usage

The usage of Yield is demonstrated in the crop-aphid-fungus model.