Julia for Biomathematics

Artwork by Luxor

Overview

Scientific computing software plays a major role in biomathematics for many reasons some of which include:

The complexity of biological systems leads to models or equations that are difficult to work with analytically but can be handled reasonably well by numerical methods.
The need to synthesize theoretical models and data from experiments.
Plots and other visualizations are important for understanding and communicating scientific ideas and results.
Computing allows us to use theoretical models to conduct experiments in silico that would be difficult or impossible to conduct otherwise.

Thus, this course will make use of scientific computing software. Three commonly used free and open-source languages for scientific computing in biomathematics are R, Julia, and Python. In this course, we will mostly use R, but it’s worth knowing at least a little about scientific computing in all three languages as each of them has its own strengths and weaknesses in the context of scientific computing for biomathematics. The goal is not to teach you to become an expert in any of these languages, but rather to become familiar with how to use them to solve some basic problems in biomathematics. Thus, this document is more of a reference than a tutorial.

This page will focus on Julia. The other two languages have their own pages:

Each of these pages will have a similar structure to this one. We do not provide details on or even an introduction to the basics of programming or the structure of any of the languages. There are many resources available for learning these things. Instead, we focus on providing some concise examples of the use of these languages for scientific computing in biomathematics. The hope is that the reader will be able to copy and modify these examples to solve their own problems.

The section titles listed in the table of contents for the page should indicate the topics or types of problems we provide examples for. For some of the more involved code, we provide links to additional webpages that provide more details. In some cases, we provide links to potentially helpful videos or web sites where the reader can learn more.

Julia

The Julia language is at the forefront of high-performance numerical and scientific computing (Bezanson et al. 2017). You can learn more about the history of Julia here.

Functions and Packages

As in R, most of the algorithms and methods we will want to use in solving problems in biomathematics are implemented as functions in Julia. Base Julia contains a number of functions including many that are useful for scientific computing. For example, base Julia contains a function sin that implements the mathematical function $\sin (x)$ . It is called as follows:

sin(pi/4.0)

0.7071067811865475

However, most of the functions we will use are not in base Julia but rather are contained in packages. A package is a collection of functions and other objects that can be loaded into Julia and used. Julia comes with a built-in package manager that allows you to install and manage packages. You can learn about the package manager through either the package manager documentation of the following video:

Managing Packages in Julia YouTube video, watch video on YouTube.

Envinronment Setup

One thing that is a little different working in Julia (and also Python) compared to R is that it’s common practice to utilize a separate environment for each project. This is done to avoid conflicts between packages and package versions that might be used in different projects. Basically, a project consists of a folder that contains files named Project.toml and Manifest.toml. You can learn the details about Julia environments by reading the environment documentation or get a short tutorial by watching the following video:

Julia environments YouTube video, watch video on YouTube.

The DrWatson.jl package provides a convenient way to manage projects and environments in Julia. You can learn more about DrWatson.jl by reading the documentation or watching the following video:

The following code uses DrWatson.jl to activate the environment for this project which is named biomath_2024:

Show the code

# import DrWatson.jl package
using DrWatson

# activate the biomath_2024 environment
@quickactivate "biomath_2024"

[ Info: Precompiling DrWatson [634d3b9d-ee7a-5ddf-bec9-22491ea816e1]

Now that we have activated our environment, we can use the package manager Pkg.jl to see which packages are installed in the environment and their versions:

Show the code

# import Pkg.jl package
import Pkg
# list packages installed in environment
Pkg.status()

Status `~/Documents/GitHub/biomath_2024/Project.toml`
⌃ [13f3f980] CairoMakie v0.11.4
  [a93c6f00] DataFrames v1.6.1
⌃ [0c46a032] DifferentialEquations v7.12.0
⌃ [31c24e10] Distributions v0.25.104
⌅ [5b8099bc] DomainSets v0.6.7
⌃ [634d3b9d] DrWatson v2.13.0
⌃ [61744808] DynamicalSystems v3.2.3
  [b964fa9f] LaTeXStrings v1.3.1
⌃ [94925ecb] MethodOfLines v0.10.4
⌅ [961ee093] ModelingToolkit v8.75.0
⌅ [8913a72c] NonlinearSolve v3.4.0
⌃ [429524aa] Optim v1.7.8
⌃ [7f7a1694] Optimization v3.20.2
⌅ [1dea7af3] OrdinaryDiffEq v6.66.0
⌅ [91a5bcdd] Plots v1.39.0
⌃ [f2b01f46] Roots v2.0.22
⌅ [0c5d862f] Symbolics v5.14.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

The nice thing about working with Julia environments is that we can easily share the environment with others. The Project.toml and Manifest.toml files contain all the information needed to recreate the environment. If you download the project folder for biomath_2024 from the GitHub repository, then all you have to do is:

Start Julia in the biomath_2024 folder which contains the Project.toml and Manifest.toml files.
Activate the environment using DrWatson.jl as shown above.
Run the instatiate command in the package manager to install all the packages listed in the Project.toml file¹.

Some of the packages most relevant for solving problems in biomathematics include:

Plots.jl for plotting.
DifferentialEquations.jl for numerical computing with differential equations.
MethodOfLines.jl for numerical computing with partial differential equations.

Plotting

using Plots, LaTeXStrings

f(x) = x^2
plot(f,-2.0,2.0, label=L"$f(x)=x^2$")

Figure 1: Plot of

f (x) = x^{2}

g(x) = x^3
plot!(g,-2.0,2.0, label=L"$g(x)=x^3$")

Figure 2: Add plot of

g (x) = x^{3}

to plot of

f (x) = x^{2}

plot([f, g],-2.0,2.0, lw=2, lc=["steelblue" "purple"],label=[L"$f(x)=x^2$" L"$g(x)=x^3$"])

Figure 3: Plots of

f (x) = x^{2}

and

g (x) = x^{3}

on same axes with some customizations.

p_1 = plot(f,-2.0,2.0, lw=2, lc="steelblue", title=L"$f(x)=x^2$")
p_2 = plot(g,-2.0,2.0, lw=2, lc="purple", title=L"$g(x)=x^3$")

plot(p_1,p_2, layout=(1,2), legend=false)

Figure 4: Plots of

f (x) = x^{2}

and

g (x) = x^{3}

Plotting in Julia YouTube video, watch video on YouTube.

Matrices and Linear Systems

Some methods for linear algebra in Julia are provided by the LinearAlgebra.jl package². The LinearAlgebra.jl package is part of the Julia standard library so it doesn’t need to be installed but must be imported. We can import the package as follows:

using LinearAlgebra

Here’s a vector:

my_vect = [1.0,-1.0,1.0]

3-element Vector{Float64}:
  1.0
 -1.0
  1.0

Here’s a matrix:

my_matrix = [1.0 0.0 -1.0;
             3.0 1.0 2.0; 
             -1.0 1.0 2.0]

3×3 Matrix{Float64}:
  1.0  0.0  -1.0
  3.0  1.0   2.0
 -1.0  1.0   2.0

Here’s the matrix-vector product:

b_vect = my_matrix * my_vect
b_vect

3-element Vector{Float64}:
 0.0
 4.0
 0.0

Here’s the determinant of the matrix:

det(my_matrix)

-3.9999999999999996

Here are the eigenvalues and eigenvectors of the matrix:

eigen(my_matrix)

Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
 -0.7320508075688776
  2.0
  2.732050807568877
vectors:
3×3 Matrix{Float64}:
  0.236174  -0.57735  -0.495572
 -0.881412  -0.57735   0.132788
  0.409065   0.57735   0.858356

Here’s the solution to the linear system:

my_matrix \ b_vect

3-element Vector{Float64}:
  1.0
 -1.0000000000000004
  1.0000000000000002

To learn more, check out

Basic Linear Algebra in Julia YouTube video, watch video on YouTube.

Root Finding

Root finding is the process of finding one or more values $x$ such that $f (x) = 0$ ³. A common application of root finding in biomathematics is to find the steady states of a dynamical system.

There are a number of different packages that facilitate root finding in Julia. Some notable one include:

We illustrate using the Roots.jl to find all of the roots of the function $f (x) = x^{2} - 1$ on the interval $[- 2, 2]$ :

using Roots

f(x) = x^2 - 1

f_roots = find_zeros(f,-2.0,2.0)

[ Info: Precompiling Roots [f2b01f46-fcfa-551c-844a-d8ac1e96c665]

2-element Vector{Float64}:
 -1.0
  1.0

Let’s make a plot of the function and its roots:

plot(f,-2.0,2.0, lc="steelblue",lw=2, label=L"$f(x)=x^2-1$")

scatter!(f_roots,zeros(length(f_roots)), color="purple",label="roots")

Figure 5: Roots of the function

f (x) = x^{2} - 1

on the interval

[- 2, 2]

The Roots.jl package focuses on scalar functions. The following link illustrates the use of NonlinearSolve.jl package to find the roots of the system of equations $f (x, y) = (x^{2} + y^{2} - 1, x^{2} - y^{2} + 0.5)$ :

Solving a nonlinear system - this page shows an example of solving a nonlinear system in Julia using the NonlinearSolve.jl package. The example also shows off some of the features of the ModelingToolkit.jl package functionality.

Optimization

Optimization is the process of finding the extrema, that is, maximum or minimum value of a function. A common application of optimization in biomathematics is to estimate parameters for a mathematical model to fit data. One of the major challenges in such optimization problems is the potential for multiple local extrema. Here, introduce the use of basic optimization methods implemented in Julia. Additional details, examples, and applications will be covered in course lectures.

There are a number of Julia packages for optimization. Here, we focus on the use of the Optim.jl package.Here’s an example of using the Optim.jl package to find the minimum of the function $f (x) = x^{2} - 1$ on the interval $[- 2, 2]$ :

using Optim

f(x) = x^2 - 1

result = optimize(f,-2.0,2.0)

x_m = Optim.minimizer(result)

y_m = Optim.minimum(result)

plot(f,-2.0,2.0, lc="steelblue",lw=2, label=L"$f(x)=x^2-1$")

scatter!([x_m],[y_m], color="purple",label="minimum")

[ Info: Precompiling Optim [429524aa-4258-5aef-a3af-852621145aeb]
[ Info: Precompiling ArrayInterfaceGPUArraysCoreExt [f5566fad-042a-5915-b3b1-71de67883923]
[ Info: Precompiling FiniteDiffStaticArraysExt [75e56524-3a34-51de-85ea-03aa6eac4b64]
[ Info: Precompiling RootsForwardDiffExt [a63bd285-541a-50a0-818b-605bb84ce7b2]
[ Info: Precompiling ForwardDiffStaticArraysExt [b74fd6d0-9da7-541f-a07d-1b6af30a262f]

Figure 6: Minimum of the function

f (x) = x^{2} - 1

on the interval

[- 2, 2]

Here’s an example of using the Optim.jl to find the minimum of the function $f (x, y) = x^{2} + y^{2} - 1$ :

f(x) = x[1]^2 + x[2]^2 - 1

result = optimize(f,[-1.0,1.0], LBFGS(); autodiff = :forward)

Optim.minimizer(result)

2-element Vector{Float64}:
 0.0
 0.0

Differential Equations and Dynamical Systems

In this section, we provide an overview and links to detailed examples of the numerical solution of initial value problems (IVPs) for ordinary differential equations (ODEs). Additionally, the linked examples illustrate computations and visualizations of the phase plane and phase portraits for systems of ODEs. For the numerical solution of ODEs, we will use the DifferentialEquations.jl package. We will also illustrate the use of the ModelingToolkit.jl package to facilitate the specification of ODEs.

Briefly, most implementations of numerical methods for IVPs for ODEs require the user to define a function that returns the value of the derivative of the solution at a given time, that is, the right hand side of the ODEs. Further, the user must specify an initial condition, a time interval over which to solve the ODEs, and the values of any parameters.

View webpage with basic examples
Intro to Differential Equations in Julia YouTube video, watch video on YouTube.

Partial Differential Equations

Partial differential equations (PDEs) generalize ordinary differential equations (ODEs) to allow for solutions that are functions of two or more variables. Typically, in the context of biomathematics, the variables represent time and spatial coordinates. Thus, PDEs allow us to model how a biological systems change in both time and space. The examples we will see in this course will mostly be convection-diffusion type equations or systems of such equations. These PDEs have the form

$\frac{\partial c}{\partial t} = \nabla \cdot (D \nabla c - v c) + R$

Special cases include diffusion equations of which the heat equation

$\frac{\partial c}{\partial t} = D \nabla^{2} c$

is a special case. In one dimension, the heat equation is

$\frac{\partial c}{\partial t} = D \frac{\partial^{2} c}{\partial x^{2}}$

Often, in order to have a complete model, we need to specify both initial conditions and boundary conditions for the PDE problem. For example in the case of the one-dimensional heat equation, we might specify the initial temperature distribution $c (x, 0)$ and the values of the flux at the boundaries, that is $\frac{\partial c}{\partial x} (0, t)$ and $\frac{\partial c}{\partial x} (L, t)$ for $t > 0$ .

There are many approaches to the numerical solution of PDEs. We will focus on the method of lines which discretizes the spatial domain and then uses ODE solvers such as those we previously described to solve the resulting system of ODEs.

For Julia, the MethodOfLines.jl package is available to facilitate the implementation of the method of lines for PDEs. The following links illustrate the use of this package:

Numerical solution of the heat equation in Julia - this page demonstrates using the MethodOfLines.jl package to obtain and visualize numerical solutions for the heat equation in one dimension.

Random Variables

While MATH 463 primarily focuses on deterministic models, stochastic models are also important in biomathematics, see for example (Allen 2011). We won’t cover stochastic models in the course but will have a few occasions where it will be useful to add some stochasticity to a deterministic model. Many different probability distributions are available in Julia through the Distributions.jl package. In this section, we introduce some functions associated with distributions and random sampling. Primarily, we focus on sampling from a distribution and accessing the probability density or probability mass function for a distribution. Further, we will demonstrate the main functions for the normal distribution and the binomial distribution. Working with other distributions is handled similarly. Additional details, examples, and applications will be covered in course lectures.

Binomial Distributions

Recall that a binomial distribution is a discrete probability distribution that models the number of successes in a sequence of $n$ independent experiments, each asking a yes–no question, and each with its own Boolean-valued outcome: success/yes/true/one (with probability $p$ ) or failure/no/false/zero (with probability $q = 1 - p$ ). A single success/failure experiment is also called a Bernoulli trial or Bernoulli experiment and a sequence of outcomes is called a Bernoulli process; for a single trial, i.e., $n = 1$ , the binomial distribution is a Bernoulli distribution. The probability of outcomes for a binomial random variable are given by the binomial probability mass function for:

$f (k, n, p) = Pr (k; n, p) = Pr (X = k) = (\binom{n}{k}) p^{k} (1 - p)^{n - k}$

where $k$ is the number of successes, $n$ is the number of trials, and $p$ is the probability of success in each trial.

Let’s work with a binomial distribution in Julia. First, we import the Random.jl, Statistics.jl, and Distributions.jl packages in order to access relevant functions:

using Random, Statistics, Distributions

Next, we create a binomial random variable with $n = 10$ and $p = 0.5$ :

n = 10
p = 0.5
binom_rv = Binomial(n,p)

Binomial{Float64}(n=10, p=0.5)

We can now access the mean, variance, and standard deviation of the random variable:

mean(binom_rv)

5.0

var(binom_rv)

2.5

std(binom_rv)

1.5811388300841898

Notice that the mean is $n p$ and the variance is $n p (1 - p)$ , as expected. Further, the standard deviation is the square root of the variance:

sqrt(2.5)

1.5811388300841898

We can work with the probability mass function as follows:

pdf(binom_rv,4)

0.20507812500000033

Figure 7 displays the probability mass function for a binomial random variable with $n = 10$ and $p = 0.5$ .

# Generate x values (possible number of successes)
x_values = 0:n

# Calculate the corresponding PDF values
pmf_values = pdf(binom_rv, x_values)

# Plot the probability mass function
bar(x_values, pmf_values, color="steelblue", xlabel="Number of Successes", ylabel="Probability", label="Binomial(n=10, p=0.5)", legend=:topleft)

# Customize the plot
title!("Binomial Distribution PMF")

Figure 7: Probability mass function for a binomial random variable with

n = 10

and

p = 0.5

We can also sample from the binomial distribution:

rand(binom_rv,10)

10-element Vector{Int64}:
 5
 4
 7
 4
 3
 4
 5
 5
 6
 5

Let’s generate a large number of samples and plot a histogram of the results, see Figure 8.

# Generate 1000 samples
samples = rand(binom_rv,1000)

# Plot a histogram of the samples
histogram(samples, bins = -0.5:1:10.5, color="steelblue", xlabel="Number of Successes", ylabel="Count", label="Binomial(n=10, p=0.5)", legend=:topleft)

# Customize the plot
title!("Binomial Distribution Samples")

Figure 8: Histogram of 1000 samples from a binomial random variable with

n = 10

and

p = 0.5

Normal Distributions

Recall that a normal distribution is a continuous probability distribution that is symmetric about the mean, showing that data near the mean are more frequent in occurrence than data far from the mean. In graph form, normal distribution will appear as a bell curve. Normal distributions are extremely important in statistics and are often used in the natural and social sciences to represent real-valued random variables whose distributions are not known. A random variable with a Gaussian distribution is said to be normally distributed and is called a normal deviate. The probability that a normal random variable with mean $μ$ and variance $σ^{2}$ has a value between $a$ and $b$ is given by the definite integral of a normal density function defined by:

$f (x ∣ μ, σ^{2}) = \frac{1}{\sqrt{2 π σ^{2}}} \exp (- \frac{(x - μ)^{2}}{2 σ^{2}})$

The parameter $σ$ is called the standard deviation.

Let’s work with a normal distribution in Julia. First, we create a normal random variable with $μ = 5$ and $σ = 1.5$ :

mu = 5.0
sigma = 1.5
norm_rv = Normal(mu,sigma)

Normal{Float64}(μ=5.0, σ=1.5)

Figure 9 displays the probability density function for a normal random variable with $μ = 5$ and $σ = 1.5$ .

f(x) = pdf(norm_rv,x);

plot(f,0.0,10.0, label=L"\mu=5, \sigma=1.5",lw=2, color="steelblue", xlabel="x", ylabel="Density", legend=:topleft)

Figure 9: Probability density function for a normal random variable with

μ = 5

and

σ = 1.5

We can sample from a normal distribution with $μ = 5$ and $σ = 1.5$ as follows:

rand(norm_rv,10)

10-element Vector{Float64}:
 5.237940841386095
 2.7008877825536337
 4.781995158534857
 4.85011463534766
 5.492672474498308
 5.149728856835479
 5.123035228301989
 4.414436122197267
 6.396973397775456
 4.462739383928602

Let’s generate a large number of samples and plot a histogram of the results, see Figure 10.

# Generate 1000 samples
samples = rand(norm_rv,1000)

# Plot a histogram of the samples
histogram(samples, bins = 0.0:0.25:10.0,normalize=:pdf, xlabel="x", ylabel="Density", label="", legend=:topleft, color="lightblue")

# Overlay the PDF
plot!(f,0.0,10.0, label=L"\mu=5, \sigma=1.5",lw=2, color="steelblue")

Figure 10: Histogram of 1000 samples from a normal random variable with

μ = 5

and

σ = 1.5

For more on probability and statistics using Julia, we recommend Statistics with Julia and Introduction to Probability for Data Science.

Symbolic Math

Symbolic Computation in Julia YouTube video, watch video on YouTube.

References

Allen, Linda J. S. 2011. An Introduction to Stochastic Processes with Applications to Biology. Second. CRC Press, Boca Raton, FL. https://mathscinet.ams.org/mathscinet-getitem?mr=2560499.

Bezanson, Jeff, Alan Edelman, Stefan Karpinski, and Viral B Shah. 2017. “Julia: A Fresh Approach to Numerical Computing.” SIAM Review 59 (1): 65–98. https://doi.org/10.1137/141000671.

Expand for Session Info

Show the code

VERSION

v"1.9.4"

Show the code

# import Pkg.jl package
import Pkg
# list packages installed in environment
Pkg.status()

Status `~/Documents/GitHub/biomath_2024/Project.toml`
⌃ [13f3f980] CairoMakie v0.11.4
  [a93c6f00] DataFrames v1.6.1
⌃ [0c46a032] DifferentialEquations v7.12.0
⌃ [31c24e10] Distributions v0.25.104
⌅ [5b8099bc] DomainSets v0.6.7
⌃ [634d3b9d] DrWatson v2.13.0
⌃ [61744808] DynamicalSystems v3.2.3
  [b964fa9f] LaTeXStrings v1.3.1
⌃ [94925ecb] MethodOfLines v0.10.4
⌅ [961ee093] ModelingToolkit v8.75.0
⌅ [8913a72c] NonlinearSolve v3.4.0
⌃ [429524aa] Optim v1.7.8
⌃ [7f7a1694] Optimization v3.20.2
⌅ [1dea7af3] OrdinaryDiffEq v6.66.0
⌅ [91a5bcdd] Plots v1.39.0
⌃ [f2b01f46] Roots v2.0.22
⌅ [0c5d862f] Symbolics v5.14.0
Info Packages marked with ⌃ and ⌅ have new versions available. Those with ⌃ may be upgradable, but those with ⌅ are restricted by compatibility constraints from upgrading. To see why use `status --outdated`

Footnotes

Note that you have to activate the environment any time you want to work in the project. However, you only need to instatiate the project the first time you use it in order to install all the packages. After that, you can just activate the environment and start working.↩︎
Another package that provides an interesting approach to solving linear systems in Julia is LinearSolve.jl.↩︎
Note that we allow for the function $f$ to be vector-valued, in which case we are looking for solutions to a system of nonlinear equations.↩︎

Reuse

CC BY-NC-SA 4.0