Basics of visualizing mathematical models

Author

Jeet Sukumaran

1 Generating sequences of values using the `range` function

1.1 Visualizing mathematical models require values generated under or sampled from those models

This is an example of a mathematical model relating \(y\) values to \(x\):

Model: A linear model

\(y = 6x + 4\)

To visualize this in Julia, we generate sets \(x\) and \(y\) values conforming to the above relationship and pass these to the plotting function. These values may be generated systematically—a sequence generated by an algorithm—or stochastically—sampled from a probabilistic distribution. Here, we will discuss the first case.

In all cases, we have to ensure that the type and range of values conform to the assumptions of the model. In this example, The \(x\) and \(y\) values will be assumed to be, each, single-dimensional real values, \(x, y \in \mathrm{R}^1\), taking values in the entire real line. The \(y\) values here are deterministic functions of the \(x\) values, and so they can be calculated from the \(x\) values, once these have been generated.

1.2 The `range` function generates sequences of values

Range objects are algorithmic representation of ranges

The actual value returned by the range function are not actually arrays of values, but rather a special “lazy” range object, with “lazy” here being a program design term meaning that the resources (particularly memory) will be allocated only if and when needed. What this means is that you can construct very large collections of values, with large ranges and high resolution, but not pay the resource cost of holding all these values in memory at the same time.

julia> r = range(1e-10, 1e10; step = 0.001)
1.0e-10:0.001:1.0e10
julia> length(r)
10000000000001
julia> r[begin]
1.0e-10
julia> r[end]
1.0e10
julia> r[begin + 2]
0.0020000001
julia> r[end - 4]
9.999999999996e9
julia> r[1]
1.0e-10
julia> r[2]
0.0010000001
julia> r[3]
0.0020000001

To actually generate the values, use the collect function, which will “unpack” the range object and return a Vector object (a single dimensional array) with the full sequence of values.

julia> r = range(0, 10; step = 2.5)
0.0:2.5:10.0

julia> collect(r)
5-element Vector{Float64}:
 0.0
 2.5
 5.0
 7.5
10.0

We can generate a set of \(x\) values over a particular range of values, continuous or otherwise, using the range function in Julia. The range function has a number of different methods, but for our purposes we will focus on specifying the range in terms of a three parameters: a start value, a stop value, and either a step size (interval between consecutive values in the sequence),

julia> range(0, 10; step = 1)
0:1:10
julia> range(0, 10; step = 0.2)
0.0:0.2:10.0
julia> range(-10, 10; step = 2)
-10:2:10

or length (number of values) in the sequence,

julia> range(0, 10; length = 11)
0.0:1.0:10.0
julia> range(0, 10; length = 5)
0.0:2.5:10.0
julia> range(-1, 1; length = 20)

The range function has a short-hand syntax, start:step:stop:

julia> 0.0:2.5:10.0
0.0:2.5:10.0

julia> collect(0.0:2.5:10.0)
5-element Vector{Float64}:
  0.0
  2.5
  5.0
  7.5
 10.0

julia> 0.0:2.5:10.0 == range(0, 10; step = 2.5)
true

1.3 Calculating and visualizing function values over a sequence

map: to apply a function to every element in a collection

The map function applies a given function to every element of a collection, returning a new collection of the results. map(x -> 6x + 4, x_values) passes each value in x_values through the function (here an anonymous function, defined using the arrow function definition syntax, x -> 6x + 4), one at a time, producing a corresponding output value for each. This is a common and idiomatic pattern in Julia: rather than writing an explicit loop, you map a transformation over a collection. The result has the same length and order as the input.

Anonymous vs. named functions

Note that here we passed map an anonymous function.

An anonymous function is a function defined without a name, typically for immediate or one-time use, using the arrow syntax, which list the input arguments of the left of the -> and output on the right. For example,x -> x^2 is an anonymous function that squares its input, while (n_apples, price_of_apple) -> n_apples * price_of_apple calculates the price of a collection of apples.

Anonymous functions are most commonly used when a function is simple and only needed in one place, defining it anonymously is cleaner than giving it a name. Typically, as functions get more complex, or for documentation and replicability reasons, we often prefer to use named functions.


model1(x) = 6x + 4

function model2(x) 
  x^2 + 2x + 3
end

map(model1, x_values)
map(model2, x_values)

We calculate the \(y\) values corresponding to each \(x\) value generated by the range object using the map function, and visualize the relationship between the \(x\) and \(y\) points using the scatter function from the CairoMakie library:

using CairoMakie
x_values = -10:0.01:10
y_values = map(x -> 6x + 4, x_values)
scatter(x_values, y_values)

or, if we do not need to store the \(y\) values beyond the visualization:

using CairoMakie
X = -10:0.01:10
scatter(X, map(x -> 6x + 4, X))

Exercise 1

Visualize larger ranges and higher resolutions (smaller steps) of \(x\)-values of this same model

\[ y = 6x + 4 \]

to see if the relationships change:

\(-1000 < x < 1000\), with step sizes of 1.0
\(0 < x < 1.0\), with step sizes of 0.001

Exercise 2

Visualize the following models in terms of the relationship between the varibles under the specified ranges:

\(f(x) = \frac{1}{x}, x \in \{1, 2, \ldots, 100\}\)
\(f(x) = x^2, -10 \le x \le 10\)

2 Plotting containers: `Figure` and `Axis` objects

Proficiencies

For pilot and exploratory work: plotting without explicit Figure and Axis objects.
Creating a Figure container object to hold all the plots.
Creating an Axis object onto which you will do the plotting, and placing it in a Figure object.
Using a stock plotting functions to plot given 2-dimensional data (\(x-\) and \(y\)-coordinates, for example).
Saving the figure to a file or displaying it in-situ.

2.1 Plotting with implict `Axis` and `Figure`


# Load the graphical module
using CairoMakie

# Generate x-values from 1 to 100
x_vals = 1:100

# Generate y-values from function applied to x-values
y_vals = map(x -> 3x + 2, x_vals)

# Plot (and display) using the `scatter` function
render = scatter(x_vals, y_vals)

# Save to file
# save("scatter1.pdf", render)
# Show figure in the REPL
# render.figure

2.2 Plotting with explicit `Axis` and `Figure` management

# Load the graphical module
using CairoMakie

# Create a `Figure` object
figure = Figure();

# Create an `Axis` object, placing it in the center of the figure.
# `[1, 1]` specifies cell (1, 1) of a grid that grows dynamically as 
# cells are specified. Here, as there is only one cell defined, it
# effectively places the `Axis` object in the center of the `Figure`
# object
ax = Axis(figure[1, 1])

# Plot some data onto the Axis object

# x-coordinate values: a sequence from 0.0 to 100.0, in 0.01 step increments
x_vals = 0:0.1:100.0

# y-coordinate values: a function of the x-coordinate values
#
# $$
# y = 4 x^2 + 2 x + 3
# $$
#
y_vals = map(x -> 4(x^2) + 2x + 3, x_vals)

# Plot using the `scatter!` function
# Note: the name is `scatter!` not `scatter`, and we pass 
# in the Axis object `ax` as the first argument
scatter!(ax, x_vals, y_vals)

# Save to file
# save("scatter1.pdf", render)
# Show figure in the REPL
# render.figure

Scatter{Tuple{Vector{Point{2, Float64}}}}

3 Visualizing the single-strain parasite virulence-transmission balance fitness landscape

Proficiencies

Generating a collection of values with range (or start:step:stop).
Computing function outputs over a collection with map.
Plotting relationships between two collections of values.
Creating an explicit Figure and Axis, and adding multiple curves to the same axis with ...! plotting functions.
Visualizing a 2-parameter model as a surface by evaluating the model on a 2D grid.

3.1 Parasite fitness function

The single-strain parasite fitness model is given as (Frank 1996):

As discussed in ?@sec-markov-property, the defining condition only depends on the present state.

\[ \begin{aligned} w = z (1 - \alpha z) \end{aligned} \tag{1}\]

where:

\(w\) is a measure of parasite fitness,
\(z\) is a measure of the parasite host exploitation severity, correlating positively with parasite reproductive and host colonization capacity, and
\(\alpha z\) is a measure of virulence, the damage a parasite does to the host as a result of exploitation \(z\), with \(\alpha\) a scaling factor, weighting the damage of the parasite reproductive activity in reducing the host resources available to the parasite.

3.2 Behavior under adaptive evolution

Let us consider host exploitation severity, \(z\), as a heritable trait that affects transmission success and virulence, and let \(w(z, \alpha)\) be the fitness of of a parasite with host exploitation severity trait value \(z\), given some virulence impact weight \(\alpha\).

If parasites vary in \(z\) and this variation is heritable, then strategies with larger \(w(z)\) tend to leave more descendants (i.e., spread more effectively).
In this simplified model, we therefore expect evolution to favor values of \(z\) that make \(w(z)\) as large as possible, so we look for the value of \(z\) that maximizes \(w(z)\).

Conceptually, we can treat \(z\) as a control or tuning variable: it is the parasite strategy we imagine changing. The parameter \(\alpha\) determines how strongly exploitation reduces host resources (how costly virulence is per unit increase in \(z\)). For a fixed \(\alpha\), plotting \(w\) against \(z\) shows a trade-off: increasing \(z\) raises transmission potential, but beyond a point the virulence cost dominates and \(w\) declines.

3.3 Visualization of fitness relationship with host exploitation severity

We will visualize the single-strain fitness model

\[ w(z;\alpha) = z(1-\alpha z), \]

where:

\(z\) is host exploitation severity (the model’s adjustable trait variable),
\(\alpha\) is a fixed scaling parameter that determines how strongly exploitation reduces transmission opportunity,
\(w\) is parasite fitness (new infected hosts per infected host under the model).

In Julia, we can generate a set of \(z\) values, evaluate \(w(z;\alpha)\), and plot \(w\) versus \(z\) for one or several choices of \(\alpha\). These values can be generated systematically, following some scheme or algorithm, or stochistically, such as by sampling from a probability distribution or model. Here, we use methods of the range function to generate the sets of values systematically at regular intervals such that they form a “grid” spanning the range of values of interest.

When we plot \(w\) versus \(z\) for a fixed \(\alpha\), we typically restrict \(z\) to

\[ 0 \le z \le \frac{1}{\alpha}, \]

so that the transmission-opportunity term \((1-\alpha z)\) stays nonnegative.

3.4 Single \(\alpha\), implicit plot (quick exploratory)

using CairoMakie

# Fitness model
w(z, α) = z * (1 - α*z)

# Choose a single α value
α = 0.5

# Generate z-values over the biologically meaningful range
z_vals = range(0.0, 1/α; length = 400)

# Evaluate the model on the sequence
w_vals = map(z -> w(z, α), z_vals)

# Quick plot with implicit Figure/Axis
lines(z_vals, w_vals)

This produces one curve showing how fitness changes as host exploitation severity \(z\) increases, holding \(\alpha\) fixed.

Multiple \(\alpha\) values on the same plot (explicit Figure and Axis) using CairoMakie

# Fitness model
w(z, α) = z * (1 - α*z)

# Choose several α values to compare
α_values = [0.2, 0.5, 1.0]

# Create Figure and Axis explicitly
fig = Figure()
ax = Axis(
    fig[1, 1],
    xlabel = "host exploitation severity (z)",
    ylabel = "parasite fitness (w)",
    title  = "Fitness curves for different α"
)

# Add one curve per α onto the same axis
for α in α_values
    z_vals = range(0.0, 1/α; length = 400)
    w_vals = map(z -> w(z, α), z_vals)
    lines!(ax, z_vals, w_vals, label = "α = $(α)")
end

# Add a legend
axislegend(ax)

fig

Here, \(\alpha\) plays the role of a parameter you “sweep” across several values, to see how the \(w\)–\(z\) relationship changes.

Surface plot: fitness as a function of \(z\) and \(\alpha\)

To visualize the full two-parameter relationship \(w(z;\alpha)\), we evaluate \(w\) on a 2D “grid” of \((z,\alpha)\) values and plot a surface.

using CairoMakie

# Fitness model
w(z, α) = z * (1 - α*z)

# Grid of α values (parameter axis)
α_vals = range(0.2, 1.0; length = 120)

# Grid of z values (trait axis)
# Use a common z-range for the surface; choose up to 1/max(α) so the full grid is comparable
z_max = 1 / maximum(α_vals)
z_vals = range(0.0, z_max; length = 200)

# Evaluate w on the (z, α) grid:
# W[i, j] = w(z_vals[i], α_vals[j])
W = [w(z, α) for z in z_vals, α in α_vals]

# Plot as a 3D surface
fig = Figure()
ax3 = Axis3(
    fig[1, 1],
    xlabel = "host exploitation severity (z)",
    ylabel = "virulence scaling (α)",
    zlabel = "parasite fitness (w)",
    title  = "Fitness surface w(z; α)"
)

surface!(ax3, z_vals, α_vals, W)

fig

This surface makes it visually explicit that \(w\) depends on both the trait value (\(z\)) and the environment/biology scaling parameter (\(\alpha\)).

Exercise 3

Single-curve practice (implicit plot)

Set \(\alpha = 0.25\) and reproduce the single-curve plot.

Increase the resolution by changing length = 400 to length = 2000.

Restrict the plotted domain to \(0 \le z \le \frac{1}{2\alpha}\) and compare the shape to the full range \(0 \le z \le \frac{1}{\alpha}\).

Exercise 4

Multiple-curve practice (explicit Axis, adding curves)

Plot curves for \(\alpha \in {0.1, 0.2, 0.4, 0.8}\) on the same axis.

Add a second axis below the first that plots the same curves using scatter! instead of lines!.

Change axis labels and title so that the plot reads well when exported.

Exercise 5

Surface practice (2D grid evaluation)

Change the \(\alpha\) grid to range(0.1, 2.0; length = 200) and regenerate the surface.

Replace the surface with a heatmap of \(w\) over \((z,\alpha)\) using heatmap!(…).

Modify the grid so that, for each \(\alpha\), you only include \(z \le 1/\alpha\) (the nonnegative-opportunity region) by setting \(w(z,\alpha) = NaN\) whenever \(z > 1/\alpha\), and observe how the plotted region changes.

References

Frank, Steven A. 1996. “Models of Parasite Virulence.” The Quarterly Review of Biology 71 (1): 37–78. https://doi.org/10.1086/419267.

1 Generating sequences of values using the range function