#head.html

#header.html#models#Cereal aphid-fungus model

Cereal aphid-fungus model

xxx

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

xxx

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.

xxx

xxx

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

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 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:

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:

xxx

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:

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:

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.

xxx

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.

xxx

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:

xxx

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.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.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.

xxx

xxx

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 = TRUE
.stopIteration = !question[accepted]
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`

has been set `TRUE`

. 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:

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`

.

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 it. 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.

Since we 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 only 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.

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=32,768\) to \(N=65,536\) for some of the nine outcome variables (arranged in rows). We judge that \(N=131,072\) is a sufficient sample size for our analysis.

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. You can tell a long story based on this sensitivity analysis. For that we refer to our paper (Saussure et al. 2023).

We already saved the data frame by setting `saveDataFrame`

to `TRUE`

in the `OutputR`

box. Since the bootstrap estimates of the confidence limits on the sensitivity indices can be somewhat time consuming to compute, those estimates are always, automatically, saved as well. This allows you to restore the analysis quickly later on.

The two files are both saved as R data frames written to your output folder. To find out where that is, use the get folders command:

```
> get folders
```

Visit that folder and you will find these two, latest files:

*biocontrol-model-sa_nnnn.Rdata**biocontrol-model-sa_nnnn-S.Rdata*

where *nnnn* is some running number. The output folder is not a safe place for files (you should empty it now and then). Copy those two files to another folder for safe-keeping.

You must write your own R script to process the two files but there is an R script to get you started residing in the same folder as the `biocontrol-model-sa.box`

script:

```
# biocontrol-model-sa-restore.R⏷
rm(list=ls(all=TRUE))
# Load standard scripts
source("<path to your input folder>/scripts/begin.R")
source("<path to your input folder>/scripts/begin-sobol.R")
# Load your data
setwd("<path to your own folder>")
load("biocontrol-model-sa_nnnn.Rdata")
load("biocontrol-model-sa_nnnn-S.Rdata")
```

The first line clears R memory just to be certain that we are working on the restored data. The scripts begin.R⏷ and begin-sobol.R⏷ are standard scripts with various useful functions found in the `scripts`

sub-folder inside your input folder. Find out where your input folder is by the get folders command:

```
> get folders
```

The two calls of `load`

will create the data frames `sim`

and `S`

in R. You can inspect them and find them pretty self-explanatory, except that in the `S`

data frame, we recommend using the bootstrap-estimated confidence limits (`LowerPercentile`

and `HigherPercentile`

) rather than the estimated standard error (`EffectSE`

).

You can treat these data however you want but the final lines in the biocontrol-model-sa-restore.R⏷ script shows you how to re-create the three plots produced by your original biocontrol-model-sa.box⏷ script:

```
...
# Prepare analysis
sobol_k = length(unique(S$Input)) - 1
sobol_N = nrow(sim)/(sobol_k+2)
iterationColumn = "iteration"
stepColumn = "step"
# Plot trends
plot_facetted(
df = sim[1:min(nrow(sim),1000),],
x = input_names(),
y = output_names(),
ncol = 1,
dir = 'h') + geom_smooth(colour='yellow')
# Plot convergence
plot_sobol_convergence()
# Plot indices
ggarrange(
plotlist=plot_sobol_indices(S)
)
```

xxx

xxx

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 of either your
`sim`

or`S`

data frame - 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 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.

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.

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

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

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}\).

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).

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

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

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

You create this figure by running the boxscript,

```
> run models/aphid/figure-10.box
```

to get the result:

You create this figure by running the boxscript,

```
> run models/aphid/figure-11.box
```

to get the result:

xxx

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.

xxx

xxx

#plugins/aphid/aphidimmigration.html

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 |

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:

xxx

xxx

#plugins/aphid/aphidindex.html

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\) |

The `AphidIndex`

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

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.

The usage of `AphidIndex`

is demonstrated in the crop-aphid-fungus model.

xxx

xxx

#plugins/aphid/aphidjuvenilesurvival.html

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 |

The `AphidJuvenileSurvival`

model computes the daily survival rate of nymphs from the inputs `cropGrowthStage`

and `temperature`

:

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

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:

xxx

xxx

#plugins/aphid/AphidNetReproduction.html

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})\) |

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} $$The usage of `AphidNetReproduction`

is demonstrated in the crop-aphid model.

xxx

xxx

#plugins/aphid/aphidoffspring.html

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] |

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:

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

The usage of `AphidOffspring`

is demonstrated in the crop-aphid model.

xxx

xxx

#plugins/aphid/biocontrol.html

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 |

The `Biocontrol`

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

The usage of `Biocontrol`

is demonstrated in the biocontrol model.

xxx

xxx

#plugins/aphid/cadaverconversion.html

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 |

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.

The usage of `CadaverConversion`

is demonstrated in the crop-aphid-fungus model.

xxx

xxx

#plugins/aphid/cadavertime.html

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 |

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:

The usage of `CadaverTime`

is demonstrated in the crop-aphid-fungus model.

xxx

xxx

#plugins/aphid/cropgrowthstage.html

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)\) |

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

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

The usage of `CropGrowthStage`

is demonstrated together with AphidImmigration.

xxx

xxx

#plugins/aphid/cropisgrowing.html

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? |

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.

The usage of `CropIsGrowing`

is demonstrated in the crop model.

xxx

xxx

#plugins/aphid/infectionrate.html

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] |

Infection spreads at the instantaneous rate

The usage of `InfectionRate`

is demonstrated in the crop-aphid-fungus model.

xxx

xxx

#plugins/aphid/prevalence.html

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 |

The `Prevalence`

box calculates the current percentages of `exposed`

aphids and of `cadavers`

.

The usage of `Prevalence`

is demonstrated in the crop-aphid-fungus model.

xxx

xxx

#plugins/aphid/yield.html

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] |

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):

The usage of `Yield`

is demonstrated in the crop-aphid-fungus model.

#right.html

Download the latest version with the newly published Cereal Aphid-Fungus model. Also includes the Virtual Greenhouse model.

5 Oct 2023

Read our paper on the Cereal Aphid-Fungus model and study the detailed documentation. Any questions? Write us.

2 Aug 2023

We remain candy-coloured until further notice.

1 Aug 2023

Any questions concerning our models and tools? Interested in visiting the lab? Want to chat online? Write us.

#footer.html