---
title: "Network Science: From Theory to Biological Applications"
subtitle: "BC204 Hands-on Workshop · Spring 2026"
author: "Savvas Paragkamian"
date: today
format:
html:
toc: true
toc-float: true
toc-depth: 3
code-fold: show
code-tools: true
theme: flatly
number-sections: true
fig-width: 8
fig-height: 5
embed-resources: true
execute:
warning: false
message: false
bibliography: ../Bibliography/Network Science Workshop 2018.bib
---
```{r}
#| label: setup
#| code-fold: true
required <- c("igraph", "poweRlaw", "ggplot2", "dplyr", "tidyr", "tibble", "purrr")
missing <- required[!required %in% installed.packages()[, "Package"]]
if (length(missing)) install.packages(missing, repos = "https://cloud.r-project.org")
library(igraph)
library(poweRlaw)
library(ggplot2)
library(dplyr)
library(tidyr)
library(tibble)
library(purrr)
library(knitr)
theme_net <- function() {
theme_bw() +
theme(panel.grid.minor = element_blank(),
strip.background = element_blank())
}
```
---
| Time | Part | Topic |
|------|------|-------|
| 0:00 – 0:20 | Part 0 | Introduction, Data structures, network types, complexity theory |
| 0:20 – 0:55 | Part 1 | Graph models (ER, BA, SW) |
| 0:55 – 1:35 | Part 2 | *C. elegans* connectome, centralities, power-law |
| 1:35 – 2:20 | — | **Break (45 min)** |
| 2:20 – 2:50 | Part 3 | Community detection (Leiden), assortativity |
| 2:50 – 3:30 | Part 4 | *E. coli* PPI, robustness, network inference |
| 3:30 – 4:00 | Part 5 | Networks as computation, GNNs |
---
# Part 0 · Data Structures and Network Types `[~20 min]`
> **Goal**: Build the vocabulary and the R objects we will use for the rest of the day. Every network analysis starts here.
Also each one kindly share:
* Name
* Expertise background
* Favourite bioinformatics tool
* AI toolkit
## Complicated ≠ Complex
These words are used interchangeably but they are not the same, and the distinction determines what kind of model is appropriate.
A **complicated** system has many parts, but the parts compose trivially. A mechanical watch has hundreds of cogs and springs, each with a specific job. Take it apart, document every piece, and you have a complete description. Behaviour is the sum of parts. Complicated systems need detailed blueprints.
A **complex** system has parts that *interact*, with feedbacks and nonlinearities and emergent properties. The interesting behaviour lives in the couplings, not the components. A three-species food web can show oscillations and regime shifts that no single species possesses. Complex systems often admit strikingly simple mathematical descriptions — because the richness is generated by a few interaction rules, not by idiosyncratic detail [@Strogatz2001a; @Kauffman1993a].
Start with the simplest possible population model. **Malthusian growth** assumes a population grows proportionally to its current size — no resource limits, no interactions:
$$N_{t+1} = r \cdot N_t \qquad \Longrightarrow \qquad N_t = N_0 \cdot r^t$$
```{r}
#| label: malthusian
#| fig-cap: "Malthusian dynamics: three growth rates, exponential trajectories. No interactions — the system is complicated only if N₀ is large, but never complex."
malthus <- function(r, n0 = 1, steps = 20) {
tibble(t = 0:steps, N = n0 * r^(0:steps), r = as.character(r))
}
map_dfr(c(0.8, 1.0, 1.2), malthus) |>
ggplot(aes(x = t, y = N, colour = r)) +
geom_line(linewidth = 1) +
scale_colour_manual(values = c("0.8" = "steelblue", "1" = "grey50", "1.2" = "coral2"),
labels = c("r = 0.8 (decline)", "r = 1.0 (steady)", "r = 1.2 (growth)")) +
labs(x = "Time step", y = "Population N", colour = NULL,
title = "Malthusian (exponential) growth") +
theme_net()
```
Malthusian dynamics is **not complex** — it is entirely predictable. The behaviour of the whole is the behaviour of one individual, scaled. No emergence is possible.
Adding a single interaction term — **density dependence**, the idea that crowding limits growth — transforms the model entirely:
$$x_{t+1} = r \cdot x_t (1 - x_t)$$
The $(1 - x_t)$ term makes each individual interact with the rest of the population. One interaction rule. One parameter. This model was proposed by ecologist Pierre Verhulst who introduced the carrying capacity in the growth model.
```{r}
#| label: logistic-map
#| fig-cap: "Logistic map bifurcation diagram: stable fixed point → period doubling → chaos. The only change from the Malthusian model is the (1 − x) term."
logistic_attractor <- function(r, n_burn = 500, n_keep = 300, x0 = 0.5) {
x <- x0
for (i in seq_len(n_burn)) x <- r * x * (1 - x)
xs <- numeric(n_keep)
for (i in seq_len(n_keep)) { x <- r * x * (1 - x); xs[i] <- x }
xs
}
r_vals <- seq(2.5, 4.0, length.out = 800)
map_dfr(r_vals, function(r) tibble(r = r, x = logistic_attractor(r))) |>
ggplot(aes(x = r, y = x)) +
geom_point(size = 0.05, alpha = 0.3, colour = "steelblue") +
labs(x = "Growth rate r", y = "Long-run population x",
title = "Logistic map bifurcation diagram") +
theme_net()
```
Beyond $r \approx 3.57$ the orbit never settles. The system is **deterministic** — the same $r$ and $x_0$ always give the same trajectory — yet **unpredictable**: two trajectories starting at almost the same $x_0$ rapidly diverge.
```{r}
#| label: chaos-sensitivity
#| fig-cap: "Sensitivity to initial conditions at r = 3.9: two trajectories differing by 0.0001 are indistinguishable for ~15 steps, then diverge completely. This is the hallmark of deterministic chaos."
logistic_ts <- function(r, x0, steps = 60) {
x <- numeric(steps + 1); x[1] <- x0
for (i in seq_len(steps)) x[i + 1] <- r * x[i] * (1 - x[i])
tibble(t = 0:steps, x = x)
}
r_chaos <- 3.9
bind_rows(
logistic_ts(r_chaos, x0 = 0.5000) |> mutate(start = "x₀ = 0.5000"),
logistic_ts(r_chaos, x0 = 0.5001) |> mutate(start = "x₀ = 0.5001")
) |>
ggplot(aes(x = t, y = x, colour = start)) +
geom_line(linewidth = 0.8) +
scale_colour_manual(values = c("steelblue", "coral2")) +
labs(x = "Time step", y = "Population x", colour = "Initial condition",
title = "Deterministic chaos: two nearly identical initial conditions",
subtitle = paste0("r = ", r_chaos, " | Δx₀ = 0.0001")) +
theme_net()
```
Beyond ~3.57 the orbit never settles — deterministic chaos from a deterministic rule. **This is emergence**: qualitatively new behaviour that is not present in any single part.
## Reaction-diffusion and spatial pattern formation
The logistic map showed how a single interaction rule generates temporal complexity. Turing's 1952 paper *The Chemical Basis of Morphogenesis* [@Turing1952] showed the same principle operating in **space**: two interacting morphogens with different diffusion rates can spontaneously break a spatially uniform state into a periodic pattern — stripes, spots, labyrinths — from nothing but random noise.
The mechanism is **short-range activation, long-range inhibition**. The activator promotes itself but also produces an inhibitor. The inhibitor diffuses faster, suppressing the activator at distance. The result is a spatial competition whose outcome depends entirely on the ratio of diffusion rates, not on any pre-existing spatial information.
This is modelled by a pair of coupled PDEs. The Gray-Scott formulation [@Pearson1993] uses two chemical species $u$ (substrate, consumed) and $v$ (catalyst, self-replicating):
$$\frac{\partial u}{\partial t} = D_u \nabla^2 u - uv^2 + F(1-u)$$
$$\frac{\partial v}{\partial t} = D_v \nabla^2 v + uv^2 - (F+k)v$$
```{r}
#| label: turing-sim
#| fig-cap: "Gray-Scott reaction-diffusion on a 80×80 grid after 5000 steps. Spots emerge from a uniform initial state perturbed only by noise — spatial complexity from two interaction rules."
#| cache: true
set.seed(42)
n <- 80
Du <- 0.16; Dv <- 0.08; F_rate <- 0.035; k <- 0.060
u <- matrix(1, n, n)
v <- matrix(0, n, n)
# Seed a central patch with inhibitor + noise
r <- 8L
cx <- (n %/% 2 - r):(n %/% 2 + r)
u[cx, cx] <- 0.50 + matrix(runif((2*r+1)^2, -0.05, 0.05), 2*r+1)
v[cx, cx] <- 0.25 + matrix(runif((2*r+1)^2, -0.05, 0.05), 2*r+1)
# Discrete Laplacian with periodic boundary conditions
lap <- function(m) {
rbind(m[-1L, ], m[1L, ]) + # up
rbind(m[n, ], m[-n, ]) + # down
cbind(m[, -1L], m[, 1L]) + # left
cbind(m[, n ], m[, -n ]) - # right
4 * m
}
dt <- 1.0
for (i in seq_len(5000L)) {
uvv <- u * v * v
u <- u + dt * (Du * lap(u) - uvv + F_rate * (1 - u))
v <- v + dt * (Dv * lap(v) + uvv - (F_rate + k) * v)
}
expand.grid(x = seq_len(n), y = seq_len(n)) |>
mutate(u = as.vector(u)) |>
ggplot(aes(x = x, y = y, fill = u)) +
geom_tile() +
scale_fill_distiller(palette = "RdYlBu", name = "u") +
coord_fixed() +
labs(title = "Turing pattern — Gray-Scott reaction-diffusion",
subtitle = paste0("Du=", Du, " Dv=", Dv,
" F=", F_rate, " k=", k)) +
theme_void() +
theme(legend.position = "right")
```
The pattern is not encoded anywhere in the initial conditions or the parameters individually. It is an **emergent property of the interaction structure** — the coupling between the two equations. Change the ratio $D_v/D_u$ and the pattern shifts from spots to stripes to labyrinths. Change it enough and the pattern disappears entirely. This is the same logic as the logistic map, now distributed across space.
A **cell** is both complicated *and* complex. Billions of molecular machines, each with specific evolved structure (complicated) interacting through feedback loops and signalling cascades that produce homeostasis, differentiation, and cell fate decisions (complex). This is why whole-cell modelling is hard in a way that aircraft engineering and magnetism are not: we cannot get away with just the blueprint, and we cannot get away with just the coarse-grained dynamics.
| System | Complicated | Complex | Implication |
|---|:---:|:---:|---|
| Mechanical watch | ✓ | — | Detailed blueprint is the right model |
| Logistic map / Ising model | — | ✓ | Simple equation captures everything |
| Epidemic, food web | — | ✓ | A few interaction rules suffice |
| Cell, brain, ecosystem | ✓ | ✓ | Neither alone works — genuinely hard |
**Networks are the natural language of complex systems.** They make interactions explicit: nodes are the parts, edges are the couplings. The same formalism applies to neurons, proteins, species, and — as we will see in Part 5 — artificial intelligence.
## The cell is not a bag of genes
JBS Haldane described the cell as a "bag of genes". We know differently now. The first things we encounter inside a cell are not genes but **molecular machines** — elaborate, evolved structures that do specific physical work: the ribosome translates mRNA, the spliceosome edits it, myosin pulls on actin, the nuclear pore complex controls everything entering or leaving the nucleus [@Hartwell1999a; @Alon2003].
Some machines are well-modelled. Ribosome queueing along mRNA has been studied since the 1960s using the **TASEP** (Totally Asymmetric Simple Exclusion Process) — a particle-exclusion model on a graph — borrowing from statistical physics. Myosin has its own lineage stretching back to Huxley's 1957 cross-bridge theory.
Many others are strikingly under-modelled despite spectacular cryo-EM structures:
```{r}
#| label: modelling-gap
tribble(
~Machine, ~Subunits, ~`Known structure`, ~`Mathematical models (approx.)`, ~`Key function`,
"Ribosome", "~80", "Yes (Å resolution)", "Thousands (TASEP, stochastic)", "mRNA translation",
"Myosin II", "~6", "Yes", "Hundreds (cross-bridge theory)", "Muscle contraction",
"Spliceosome", ">100", "Yes (cryo-EM)", "~Tens", "pre-mRNA splicing",
"Chaperonin (GroEL)", "~14", "Yes", "~Tens", "protein folding",
"V-ATPase", "~30", "Yes", "~Tens", "proton pumping",
"Nuclear pore complex", "~120", "Yes (cryo-ET)", "Very few", "nuclear transport"
) |> kable(caption = "Molecular machines: structural knowledge vs mathematical models (after Stumpf)")
```
The challenge is not obtaining structures — cryo-EM has transformed that. The hard step is translating **structure → dynamics → predictive mathematical model**. Fewer than 1% of published models ever reach shared repositories like BioModels.
> **Q**: The spliceosome has over 100 subunits and near-atomic resolution structures — yet almost no mechanistic models. Is this a complicated problem, a complex one, or both? What would a network model of the spliceosome represent?
## The human is not determined by gene count
When the human genome was sequenced, the expectation was ~100,000 genes — one for each protein, each trait, each behaviour. The actual number was ~20,000–25,000. About the same as a nematode worm.
```{r}
#| label: gene-count
#| fig-cap: "Gene count does not track organismal complexity. A rice plant has more protein-coding genes than a human. The complexity of an organism lives in its regulatory interactions, not its parts list."
tribble(
~Organism, ~`Common name`, ~`Genome size (Mb)`, ~`Protein-coding genes`,
"Homo sapiens", "Human", 3200, ~20500,
"Mus musculus", "Mouse", 2800, ~20000,
"Drosophila melanogaster","Fruit fly", 180, ~14000,
"Caenorhabditis elegans","Nematode worm", 100, ~20000,
"Arabidopsis thaliana", "Thale cress", 135, ~27000,
"Oryza sativa", "Rice", 430, ~38000,
"Paris japonica", "Canopy plant", 149000, ~30000
) |>
arrange(desc(`Protein-coding genes`)) |>
kable(caption = "Gene count does not track organismal complexity (G-value paradox)")
```
The resolution is not in the genes but in the **interactions**. Alternative splicing, post-translational modification, and regulatory networks mean that 20,000 genes can generate far more distinct states than a parts list suggests. A network of $N$ nodes with binary states has $2^N$ possible configurations. The cell navigates this state space through its interaction structure — not through gene number.
## Ecosystem function cannot be predicted from species count
The same logic extends to ecology. The number of species in a community — biodiversity — does not determine ecosystem function. A species that appears ecologically redundant under current conditions may be the only one capable of maintaining function under novel conditions (drought, invasion, temperature shift). Remove it and the system can shift abruptly to a different state.
The diversity–stability debate has run since MacArthur (1955) and May (1972). May showed analytically that **random complexity destabilises**: a randomly wired ecological network with $S$ species and connectivity $C$ is stable only if $SC\sigma^2 < 1$ (where $\sigma^2$ is interaction-strength variance). More species and more connections, holding strength constant, makes the system *less* stable. This is the opposite of the intuitive view.
The resolution is that real ecological networks are not random. They have structure: weak interactions dominate, strong interactions are rare and embedded in stabilising motifs, interaction topology shapes stability more than species richness alone [@McCann2000; @Bascompte2009a; @Pascual2007; @Montoya2006].
## More is Different
In 1972, Philip W. Anderson published a one-page paper in *Science* with the title *More is Different* [@Anderson1972]. It is one of the founding texts of complexity science. The same principle runs through Kauffman's work on self-organisation in evolution [@Kauffman1993a].
The argument is simple. The reductionist programme — understand the parts, and you understand the whole — works within a level of organisation. Particle physics is not wrong about quarks. Chemistry is not wrong about molecules. But knowing the laws of quarks does not give you the laws of chemistry. Knowing the laws of chemistry does not give you the laws of cell biology. At each level, **new properties appear that are not present at the level below and cannot be derived from it by straightforward extrapolation**.
Anderson called this *broken symmetry*: macroscopic systems do not share the symmetries of their microscopic constituents. A crystal has a preferred orientation; the underlying quantum mechanics does not. A flock has a collective direction; individual birds do not. A gene regulatory network has bistability; individual protein interactions do not. Prigogine extended this idea into thermodynamics: far-from-equilibrium systems spontaneously break time symmetry, generating irreversible order — dissipative structures — from fluctuations [@Prigogine1984; @Antoniou1997].
The practical consequence: you cannot compress a complex system by reducing it. The emergent properties are not encoded in the parts — they are encoded in the **relational structure** between parts. This is what networks represent.
## Reductionism, holism, and where network science fits
Science has operated mostly under a reductionist assumption: find the smallest parts, characterise them precisely, and reconstruct the whole. This works exceptionally well for some problems — the periodic table, protein structure, gene function. It fails for others — why a particular ecosystem collapses, how a cell decides to differentiate, whether a drug will have off-target effects in a living network of 20,000 interacting genes.
Holism — the whole has properties not present in the parts — is often dismissed as unscientific because it seems to resist quantification. The contribution of complex systems science, and network science specifically, is to make holism **precise and computable**:
| Approach | Core assumption | Right tool | Fails when |
|---|---|---|---|
| Reductionism | Parts determine whole | Structural biology, genetics | Emergent behaviour dominates |
| Holism (naive) | Whole transcends parts | — | No predictive mechanism |
| Network science | Interactions determine behaviour | Graph theory, dynamics | Interactions unknown or unmeasurable |
Network science is not anti-reductionist. It is a framework for moving between levels: the nodes are the output of reductionist analysis (here is the protein, here is the neuron), and the edges encode the interactions that produce emergent behaviour at the system level [@Boccaletti2006; @Newman2010d; @Sole2002b; @Sole2004]. The cell is not a bag of genes. It is a network.
The formalism used in this workshop treats networks as static. Real biological networks are not: interactions are context-dependent, time-varying, and often layered across multiple types simultaneously. These extensions are formalised in temporal networks [@Holme2012] and multiplex networks [@Cozzo2017].
**Networks are the natural language of complex systems.** They make interactions explicit: nodes are the parts, edges are the couplings. The same formalism applies to neurons, proteins, species, and — as we will see in Part 5 — artificial intelligence.
## What is a graph?
A network is modelled mathematically as a **graph** $G(V, E)$: a set of $N$ nodes $V$ and a set of edges $E \subseteq V \times V$ [@Crane2018; @Wasserman1994].
Four assumptions underpin this abstraction:
1. All nodes belong to the same variable type (e.g., all proteins, or all neurons).
2. All edges belong to the same variable type (e.g., all physical interactions).
3. Each edge connects exactly two nodes — a **binary relation**.
4. Both nodes and edges are time-independent (a static snapshot).
These assumptions are often violated in real biology. Knowing *when* they break tells you when to upgrade to a richer model (temporal, multilayer, hypergraph).
## Data structures for networks
A network can be stored in three equivalent ways:
| Structure | Best for |
|---|---|
| **Adjacency matrix** | Dense networks, linear algebra, visualisation as a heatmap |
| **Edge list** | Sparse networks, tidy data pipelines |
| Incidence matrix | Bipartite networks (two node types) |
### Adjacency matrix
Rows and columns are nodes; a cell is 1 if an edge exists (0 otherwise). For a **directed** graph the matrix is asymmetric — entry $(i, j) = 1$ means an edge *from* $i$ *to* $j$. For an **undirected** graph it is symmetric.
```{r}
#| label: adjacency-matrix
set.seed(7)
n_nodes <- 5
node_names <- LETTERS[1:n_nodes]
adj <- matrix(
sample(c(0L, 1L), n_nodes^2, replace = TRUE, prob = c(0.5, 0.5)),
nrow = n_nodes, ncol = n_nodes,
dimnames = list(node_names, node_names)
)
diag(adj) <- 0L # no self-loops
kable(adj, caption = "Adjacency matrix — rows are sources, columns are targets")
```
Convert to an igraph object and plot:
```{r}
#| label: adj-to-graph
#| fig-cap: "Graph built from the adjacency matrix"
g_adj <- graph_from_adjacency_matrix(adj, mode = "directed")
plot(g_adj, vertex.color = "steelblue", vertex.label.color = "white",
edge.arrow.size = 0.6, main = "From adjacency matrix")
```
Visualise the same matrix as a heatmap — useful for spotting asymmetry (directed) or symmetry (undirected):
```{r}
#| label: adj-heatmap
#| fig-cap: "Adjacency matrix as a heatmap — asymmetry signals a directed network"
as.data.frame(adj) |>
rownames_to_column("from") |>
pivot_longer(-from, names_to = "to", values_to = "edge") |>
mutate(from = factor(from, levels = rev(sort(unique(from))))) |>
ggplot(aes(x = to, y = from, fill = factor(edge))) +
geom_tile(colour = "white") +
scale_fill_manual(values = c("0" = "grey92", "1" = "steelblue"),
labels = c("absent", "present")) +
scale_x_discrete(position = "top") +
coord_fixed() +
labs(x = "To", y = "From", fill = "Edge", title = "Adjacency matrix heatmap") +
theme_net()
```
### Edge list
An edge list stores only the existing edges — one row per edge. It is more memory-efficient than a matrix for sparse networks (most real biological networks have density < 1 %).
```{r}
#| label: edge-list
el <- as.data.frame(which(adj == 1, arr.ind = TRUE)) |>
transmute(from = node_names[row], to = node_names[col])
kable(el, caption = "Edge list derived from the same network")
```
Build a graph directly from the edge list:
```{r}
#| label: el-to-graph
#| fig-cap: "The same graph built from the edge list — identical result"
g_el <- graph_from_edgelist(as.matrix(el), directed = TRUE)
plot(g_el, vertex.color = "coral2", vertex.label.color = "white",
edge.arrow.size = 0.6, main = "From edge list")
```
### Weighted edges
Edges can carry a numeric **weight** (interaction strength, distance, correlation). Weights extend binary graphs to continuous ones.
```{r}
#| label: weighted-edgelist
el_w <- el |> mutate(weight = round(runif(n(), 0.1, 1.0), 2))
kable(el_w, caption = "Weighted edge list")
```
## Network types
Real networks differ in four properties you must check before any analysis:
| Property | Values | Example |
|---|---|---|
| **Direction** | directed / undirected | synapses (directed); co-expression (undirected) |
| **Weight** | weighted / binary | interaction strength; presence/absence |
| **Sign** | signed / unsigned | activation (+) / inhibition (−) in gene regulation |
| **Simplicity** | simple / multigraph | single edge per pair; multiple experimental evidence |
```{r}
#| label: network-type-checklist
g_w <- graph_from_data_frame(el_w, directed = TRUE)
E(g_w)$weight <- el_w$weight
tibble(
Property = c("Directed", "Weighted", "Simple", "Has loops"),
Value = c(is_directed(g_w), is_weighted(g_w),
is_simple(g_w), any(which_loop(g_w)))
) |> kable(caption = "Network type checklist for our toy graph")
```
## R objects for network data
Three representations are used across this workshop:
| Object | Package | When to use |
|---|---|---|
| `igraph` | igraph | All network analysis — fast C backend |
| `matrix` / `tibble` | base / tidyverse | Raw import, pre-processing, and statistics |
For a practical guide to network analysis in R see @Luke2015; for the statistical foundations see @Kolaczyk2009a.
```{r}
#| label: r-objects
# igraph object — the primary representation used throughout this workshop
g_ig <- graph_from_data_frame(el_w, directed = TRUE)
g_ig
```
> **Exercise 0**: Make the toy graph undirected (`as.undirected(g_ig)`). Check again: is the adjacency matrix now symmetric? Does the edge count change, and why?
---
# Part 1 · Graph Models `[~35 min]`
> **Goal**: Understand the three canonical random graph models as null hypotheses for real biological networks.
## Erdős–Rényi: the null model
The ER model places edges between *N* nodes with probability *p*, independently. It predicts a Poisson degree distribution — no hubs. Real biological networks almost never look like this, which is exactly why it is useful as a null [@Newman2010d; @Bollobas].
```{r}
#| label: er-phase-transition
#| fig-cap: "Phase transition in ER networks: a giant component suddenly appears at ⟨k⟩ = 1"
set.seed(42)
n <- 500
probs <- seq(0, 0.012, length.out = 80)
er_stats <- map_dfr(probs, function(p) {
g <- sample_gnp(n, p, directed = FALSE)
comps <- components(g)
tibble(mean_k = mean(degree(g)),
giant_frac = max(comps$csize) / n)
})
ggplot(er_stats, aes(x = mean_k, y = giant_frac)) +
geom_point(colour = "steelblue", size = 1.5) +
geom_vline(xintercept = 1, linetype = "dashed", colour = "firebrick") +
labs(x = "Mean degree ⟨k⟩", y = "Giant component / N",
title = "ER percolation threshold",
caption = "Dashed line: ⟨k⟩ = 1 — the point where a giant component emerges") +
theme_net()
```
> **Q**: A pathogen progressively disables random host proteins. Which side of this threshold matters for the host? For the pathogen's spread through a PPI network?
## Barabási–Albert: preferential attachment
Each new node attaches with probability proportional to existing degree — "rich get richer". This produces a power-law degree distribution and explains why hubs appear in PPIs, metabolic networks, and co-authorship graphs [@Barabasi1999a; @Barabasi1999c; @Albert2002; @Dorogovtsev2003].
```{r}
#| label: ba-degree
#| fig-cap: "BA model: straight line on log-log axes is the signature of a power law"
set.seed(42)
ba <- sample_pa(2000, power = 1, directed = FALSE)
tibble(degree = degree(ba)) |>
count(degree) |>
ggplot(aes(x = degree, y = n)) +
geom_point(colour = "coral2") +
scale_x_log10() + scale_y_log10() +
labs(x = "Degree (log)", y = "Count (log)",
title = "BA model degree distribution") +
theme_net()
```
## Watts–Strogatz: small-world networks
A regular lattice rewired with probability β. At intermediate β: short average path lengths (like a random graph) AND high clustering (like a lattice) — simultaneously. This is the *small-world* regime [@Watts1998; @Strogatz2001a].
```{r}
#| label: sw-regime
#| fig-cap: "Small-world window: C stays high while L drops, between β ≈ 0.01 and β ≈ 0.1"
set.seed(42)
betas <- 10^seq(-3, 0, length.out = 40)
lattice <- make_lattice(length = 200, dim = 1, nei = 4, circular = TRUE)
C0 <- transitivity(lattice, type = "average")
L0 <- mean_distance(lattice, directed = FALSE)
sw_stats <- map_dfr(betas, function(b) {
g <- sample_smallworld(1, 200, 4, b)
tibble(beta = b,
C_norm = transitivity(g, type = "average") / C0,
L_norm = mean_distance(g, directed = FALSE) / L0)
})
sw_stats |>
pivot_longer(c(C_norm, L_norm), names_to = "metric", values_to = "value") |>
ggplot(aes(x = beta, y = value, colour = metric)) +
geom_line(linewidth = 1) +
scale_x_log10() +
scale_colour_manual(values = c(C_norm = "steelblue", L_norm = "coral2"),
labels = c("C / C₀ (clustering)", "L / L₀ (path length)")) +
labs(x = "Rewiring probability β (log)", y = "Normalised value",
colour = NULL, title = "Watts–Strogatz small-world regime") +
theme_net()
```
> **Exercise 1**: Generate ER and BA networks each with N = 300 and approximately the same number of edges as a Watts–Strogatz graph with N = 300, K = 4, β = 0.05. Compare clustering coefficient and mean path length across the three. Which model is "most small-world"?
---
# Part 2 · C. elegans Connectome `[~40 min]`
The frontal neuronal network of *C. elegans*: 302 neurons (nodes), chemical synapses and gap junctions (directed edges). Fully mapped from electron microscopy — one of the few complete connectomes in biology [@Kaiser2006; @Varshney2011].
## Import and network checklist
Every network analysis should start with this checklist before computing anything.
```{r}
#| label: celegans-import
edges <- read.table("../Data/C-elegans-frontal.txt", header = TRUE, sep = " ") |>
mutate(across(c(FromNodeId, ToNodeId), as.character))
meta <- read.csv("../Data/C-elegans-frontal-meta.csv") |>
rename(x = posx, y = posy) |>
mutate(node_id = as.character(node_id))
net <- graph_from_data_frame(edges, directed = TRUE, vertices = meta)
tibble(
property = c("Nodes", "Edges", "Directed", "Simple",
"Weighted", "Connected", "Components"),
value = c(vcount(net), ecount(net), is_directed(net),
is_simple(net), is_weighted(net), is_connected(net),
components(net)$no)
) |> kable(caption = "C. elegans connectome — network checklist")
```
```{r}
#| label: giant-component
# Work with the largest weakly connected component
giant <- decompose(net, min.vertices = 10)[[1]]
cat("Giant component:", vcount(giant), "nodes |", ecount(giant), "edges\n")
```
## Spatial visualisation
Neuron positions are anatomical (inferred from EM data), so this layout is biologically meaningful — not a force-directed aesthetic choice.
```{r}
#| label: spatial-plot
#| fig-cap: "C. elegans frontal connectome in anatomical space"
layout_anat <- cbind(V(giant)$x, V(giant)$y)
plot(giant,
layout = layout_anat,
vertex.size = 3,
vertex.color = "steelblue",
vertex.label = NA,
edge.arrow.size = 0.2,
edge.color = adjustcolor("grey50", alpha.f = 0.25),
main = "C. elegans frontal connectome — anatomical layout")
```
## Degree distribution and power-law test
```{r}
#| label: degree-powerlaw
#| fig-cap: "Degree distribution: empirical (blue) vs power-law fit (red line from xmin)"
deg_all <- degree(giant, mode = "total")
deg_pos <- deg_all[deg_all > 0]
# Maximum-likelihood fit following Clauset, Shalizi & Newman [@Clauset2009]; implementation via poweRlaw [@Gillespie2014]
pl <- displ$new(deg_pos)
pl$setXmin(estimate_xmin(pl))
pl$setPars(estimate_pars(pl))
cat(sprintf("Power-law: α = %.3f | xmin = %d | tail n = %d\n",
pl$pars, pl$xmin, sum(deg_pos >= pl$xmin)))
# Empirical complementary CDF
ccdf_df <- tibble(k = sort(unique(deg_pos))) |>
mutate(p = map_dbl(k, \(x) mean(deg_pos >= x)))
# Theoretical power-law CCDF
k_range <- seq(pl$xmin, max(deg_pos))
tail_frac <- mean(deg_pos >= pl$xmin)
pl_ccdf_df <- tibble(
k = k_range,
p = tail_frac * (k_range / pl$xmin) ^ (-(pl$pars - 1))
)
ggplot(ccdf_df, aes(x = k, y = p)) +
geom_point(colour = "steelblue", size = 1.5) +
geom_line(data = pl_ccdf_df, aes(x = k, y = p),
colour = "firebrick", linewidth = 1) +
scale_x_log10() + scale_y_log10() +
labs(x = "Degree k (log)", y = "P(K ≥ k) [log]",
title = "C. elegans degree distribution",
subtitle = sprintf("Power-law fit: α = %.2f, xmin = %d", pl$pars, pl$xmin)) +
theme_net()
```
## Centralities
```{r}
#| label: centralities
node_tbl <- igraph::as_data_frame(giant, what = "vertices") |>
mutate(
degree = degree(giant, mode = "total"),
degree_in = degree(giant, mode = "in"),
degree_out = degree(giant, mode = "out"),
betweenness = betweenness(giant, directed = TRUE),
closeness = closeness(giant, mode = "all"),
pagerank = page_rank(giant)$vector
)
```
Degree measures local connectivity. Betweenness measures how often a neuron lies on the shortest path between other pairs — a bottleneck [@Freeman1979b; @Freeman1991c]. PageRank propagates influence through the directed structure. For a systematic treatment of centrality axioms see @Boldi2014a; for an overview of applications see @Lu2016a. In neural networks specifically, centrality measures have been used to identify functionally critical regions [@Bullmore2009; @Rubinov2010b].
```{r}
#| label: centrality-scatter
#| fig-cap: "Nodes above the trend are bottlenecks: high betweenness despite moderate degree"
ggplot(node_tbl, aes(x = degree, y = betweenness)) +
geom_point(colour = "steelblue", alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE, colour = "firebrick", linewidth = 0.8) +
labs(x = "Degree", y = "Betweenness centrality",
title = "Hubs vs bottlenecks in C. elegans") +
theme_net()
```
```{r}
#| label: top-neurons
node_tbl |>
arrange(desc(betweenness)) |>
select(name, degree, betweenness, closeness, pagerank) |>
head(10) |>
kable(digits = 3,
caption = "Top-10 neurons by betweenness centrality",
col.names = c("Neuron", "Degree", "Betweenness", "Closeness", "PageRank"))
```
The power-law fit should be interpreted with caution. @Stumpf2012 and @Broido2018 show that strict scale-free degree distributions are rare in empirical networks, and that random subsamples of scale-free networks are not themselves scale-free [@Stumpf2005a]. The statistical test provided by @Clauset2009 and implemented in `poweRlaw` [@Gillespie2014] is necessary but not sufficient. @Newman2005 provides a broader treatment of heavy-tailed distributions.
> **Exercise 2**: Identify neurons in the top-10 by betweenness that are *not* in the top-10 by degree. These are the biological bottlenecks. Look up one of them — what is its known function in *C. elegans* neurobiology?
---
::: {.callout-note}
## ☕ Break — 45 minutes
:::
# Part 3 · Community Structure `[~30 min]`
Neurons do not operate in isolation — they form functional modules [@Hartwell1999a; @Koch2012]. Community detection asks: which sets of neurons are more densely interconnected with each other than with the rest? For a review of algorithms and quality metrics see @Newman2006, @Newman2015, and @Emmons2016. Detected communities should not be over-interpreted: @Ingram2006a showed that network structure alone does not determine function.
## Leiden algorithm
Leiden [@Traag2019] solves a resolution-limit problem in Louvain: it guarantees that all communities are internally well-connected and produces more stable partitions. It is the current standard for biological network analysis.
```{r}
#| label: leiden-communities
# Leiden requires an undirected simple graph
g_ud <- as.undirected(giant) |> simplify()
comm_leid <- cluster_leiden(g_ud, resolution_parameter = 0.5)
V(g_ud)$community <- as.character(membership(comm_leid))
n_comm <- max(membership(comm_leid))
cat("Communities detected (resolution = 0.5):", n_comm, "\n")
```
```{r}
#| label: spatial-communities
#| fig-cap: "Community membership mapped onto anatomical positions — do communities cluster spatially?"
comm_pal <- hcl.colors(n_comm, palette = "Dark 2")
layout_anat_ud <- cbind(V(g_ud)$x, V(g_ud)$y)
plot(g_ud,
layout = layout_anat_ud,
vertex.size = 4,
vertex.color = comm_pal[as.integer(V(g_ud)$community)],
vertex.label = NA,
edge.color = adjustcolor("grey70", alpha.f = 0.12),
main = "C. elegans communities (Leiden, resolution = 0.5)")
```
```{r}
#| label: kk-communities
#| fig-cap: "Force-directed layout reveals intra-community density regardless of anatomy"
set.seed(42)
plot(g_ud,
layout = layout_with_kk(g_ud),
vertex.size = 4,
vertex.color = comm_pal[as.integer(V(g_ud)$community)],
vertex.label = NA,
edge.color = adjustcolor("grey70", alpha.f = 0.12),
main = "C. elegans communities (Kamada-Kawai layout)")
```
## Resolution sensitivity
The `resolution` parameter controls granularity: higher values split communities further. There is no universally correct value — choose based on the biological question.
```{r}
#| label: resolution-sweep
#| fig-cap: "Number of detected communities increases with resolution parameter"
resolutions <- seq(0.1, 2.5, by = 0.1)
res_df <- map_dfr(resolutions, function(r) {
n <- max(membership(cluster_leiden(g_ud, resolution_parameter = r)))
tibble(resolution = r, n_communities = n)
})
ggplot(res_df, aes(x = resolution, y = n_communities)) +
geom_line(colour = "steelblue") +
geom_point(colour = "steelblue") +
labs(x = "Resolution parameter", y = "Number of communities",
title = "Leiden resolution sensitivity") +
theme_net()
```
## Degree assortativity
Do highly connected neurons preferentially synapse onto other highly connected neurons? This tendency — assortative mixing by degree — was formalised by @Newman2010d and has been found to differ systematically between social (assortative) and biological networks (disassortative).
```{r}
#| label: assortativity
cat("Degree assortativity:", round(assortativity_degree(giant, directed = TRUE), 3), "\n")
knn_vals <- knn(giant)$knn
tibble(name = names(knn_vals), knn = knn_vals) |>
left_join(node_tbl, by = "name") |>
ggplot(aes(x = degree, y = knn)) +
geom_point(colour = "steelblue", alpha = 0.7) +
geom_smooth(method = "lm", se = FALSE, colour = "firebrick") +
labs(x = "Degree", y = "Mean neighbour degree",
title = "Degree assortativity in C. elegans") +
theme_net()
```
> **Exercise 3**: Change the Leiden `resolution` to 1.5. Do the communities still align with anatomical positions, or do they fragment? What does a negative assortativity value mean biologically in a neural circuit?
---
# Part 4 · E. coli PPI and Network Robustness `[~40 min]`
## Protein–protein interaction network
BioGRID curates experimentally validated PPIs from the literature [@Stark2006; @Krogan2006b]. We keep only *physical* interactions (direct binding or stable complex), not *genetic* interactions (functional relationships without physical contact) [@Boucher2013a; @Boone2007a]. Large-scale genetic interaction maps such as @Costanzo2016 reveal how functional buffering is distributed across the network.
```{r}
#| label: ecoli-import
biogrid_raw <- read.delim(
"../Data/BIOGRID-ORGANISM-Escherichia_coli_K12_W3110-3.5.165.mitab.txt",
sep = "\t", header = TRUE
)
physical_types <- c(
"psi-mi:MI:0407(direct interaction)",
"psi-mi:MI:0915(physical association)"
)
ppi_edges <- biogrid_raw |>
filter(Interaction.Types %in% physical_types) |>
select(from = Alt.IDs.Interactor.A, to = Alt.IDs.Interactor.B)
ppi_net <- graph_from_data_frame(ppi_edges, directed = FALSE) |>
simplify()
cat("E. coli PPI — nodes:", vcount(ppi_net),
"| edges:", ecount(ppi_net),
"| density:", round(edge_density(ppi_net), 5), "\n")
```
```{r}
#| label: ecoli-degree
#| fig-cap: "E. coli PPI degree distribution — linear on log-log, consistent with a power law"
deg_ppi <- degree(ppi_net)
tibble(degree = deg_ppi) |>
count(degree) |>
ggplot(aes(x = degree, y = n)) +
geom_point(colour = "coral2") +
scale_x_log10() + scale_y_log10() +
labs(x = "Degree (log)", y = "Count (log)",
title = "E. coli PPI degree distribution") +
theme_net()
```
## Robustness: random failures vs targeted attacks
Scale-free networks tolerate random node removal (rare hubs are unlikely to be hit) but collapse quickly under targeted hub removal [@Reka2000a; @Albert2002]. This asymmetry has direct implications for antibiotic target selection and network-based drug design. The same topological vulnerability shapes epidemic spreading: highly connected nodes act as super-spreaders [@Vespignani2012; @Barrat2008b]. For dynamical systems perspectives on network processes see @Porter2016.
```{r}
#| label: robustness
#| cache: true
#| fig-cap: "Percolation: the PPI network is robust to random failures but fragile under hub removal"
percolate <- function(g, strategy = c("random", "targeted"), n_steps = 80) {
strategy <- match.arg(strategy)
N <- vcount(g)
V(g)$id <- seq_len(N)
step_at <- round(seq(1, N - 1, length.out = n_steps))
results <- vector("list", n_steps)
for (i in seq_along(step_at)) {
comps <- components(g)
results[[i]] <- tibble(
removed_frac = (N - vcount(g)) / N,
giant_frac = max(comps$csize) / N
)
del <- if (strategy == "random") {
sample(V(g)$id, 1)
} else {
which.max(degree(g))
}
g <- delete_vertices(g, del)
V(g)$id <- seq_len(vcount(g))
}
bind_rows(results) |> mutate(strategy = strategy)
}
set.seed(42)
ppi_giant <- decompose(ppi_net)[[1]]
robust_df <- bind_rows(
percolate(ppi_giant, "random"),
percolate(ppi_giant, "targeted")
)
ggplot(robust_df, aes(x = removed_frac, y = giant_frac, colour = strategy)) +
geom_line(linewidth = 1) +
scale_colour_manual(values = c(random = "steelblue", targeted = "firebrick"),
labels = c("Random failure", "Targeted hub removal")) +
labs(x = "Fraction of proteins removed",
y = "Giant component fraction",
colour = NULL,
title = "E. coli PPI robustness") +
theme_net()
```
## Network inference: from correlation to edges
Most biological networks are not directly observed — they are *inferred* from expression data [@Ideker2012a; @Teschendorff2014; @Glass2014]. The simplest approach: correlate gene expression profiles across conditions and apply a threshold. The threshold choice is a biological decision, not a statistical one.
```{r}
#| label: inference-threshold
#| fig-cap: "Network structure is highly sensitive to the correlation threshold — there is no free lunch"
set.seed(42)
n_genes <- 60
n_samples <- 25
expr <- matrix(rnorm(n_genes * n_samples), nrow = n_genes,
dimnames = list(paste0("g", seq_len(n_genes)), NULL))
# Plant a co-expression module in genes 1–15
expr[1:15, ] <- expr[1:15, ] + rnorm(n_samples, mean = 2.5, sd = 0.5)
cor_mat <- cor(t(expr))
thresholds <- seq(0.3, 0.9, by = 0.05)
thr_stats <- map_dfr(thresholds, function(thr) {
adj <- (abs(cor_mat) > thr) * 1L
diag(adj) <- 0L
g <- graph_from_adjacency_matrix(adj, mode = "undirected")
tibble(threshold = thr,
density = edge_density(g),
n_components = components(g)$no,
mean_degree = mean(degree(g)))
})
thr_stats |>
pivot_longer(-threshold, names_to = "metric", values_to = "value") |>
ggplot(aes(x = threshold, y = value)) +
geom_line(colour = "steelblue") + geom_point(colour = "steelblue") +
facet_wrap(~ metric, scales = "free_y",
labeller = as_labeller(c(
density = "Edge density",
n_components = "Connected components",
mean_degree = "Mean degree"
))) +
labs(x = "Correlation threshold", y = NULL,
title = "Effect of threshold on inferred network structure") +
theme_net()
```
> **Exercise 4**: At what threshold does the inferred network fragment into more than 10 components? Is a denser or sparser network more biologically credible for a 25-sample experiment? What additional information (besides the threshold) would you need to validate this network?
---
# Part 5 · Networks as Computation `[~30 min]`
> **Goal**: Show that the network formalism we used all day is also the foundation of the most powerful AI tools applied to biological data. No optimisation theory or GPU systems — the bridge is relational structure and information flow.
## Artificial neural networks are graphs
An artificial neural network is a directed weighted graph. Nodes are neurons. Edges are weights. That is not a metaphor — it is the definition.
The first difference from biological networks is topological. A **multi-layer perceptron** (MLP) enforces a strict layered structure: all edges go forward, from layer $l$ to layer $l+1$, with no shortcuts, no feedback, no hubs.
```{r}
#| label: mlp-build
layer_sizes <- c(input = 4, hidden1 = 6, hidden2 = 6, output = 2)
layer_names <- names(layer_sizes)
nodes_mlp <- map_dfr(seq_along(layer_sizes), function(i) {
tibble(layer = layer_names[i],
layer_idx = i,
neuron_idx = seq_len(layer_sizes[i]),
name = paste0(layer_names[i], "_", seq_len(layer_sizes[i])))
}) |> mutate(x = layer_idx,
y = neuron_idx - (layer_sizes[layer] + 1) / 2)
edges_mlp <- map_dfr(seq_len(length(layer_sizes) - 1), function(i) {
expand.grid(from = nodes_mlp$name[nodes_mlp$layer == layer_names[i]],
to = nodes_mlp$name[nodes_mlp$layer == layer_names[i + 1]],
stringsAsFactors = FALSE)
})
g_mlp <- graph_from_data_frame(edges_mlp,
vertices = select(nodes_mlp, name, everything()), directed = TRUE)
```
```{r}
#| label: mlp-plot
#| fig-cap: "MLP as a directed graph. Every neuron in one layer connects to every neuron in the next — fully connected, no hubs, no shortcuts."
layer_pal <- hcl.colors(length(layer_sizes), palette = "Dark 2")
layer_idx <- match(V(g_mlp)$layer, layer_names)
plot(g_mlp,
layout = as.matrix(select(nodes_mlp, x, y)),
vertex.color = layer_pal[layer_idx],
vertex.size = 20,
vertex.label = V(g_mlp)$neuron_idx,
vertex.label.color = "white",
vertex.label.cex = 0.8,
edge.arrow.size = 0.3,
edge.color = adjustcolor("grey50", alpha.f = 0.3),
main = "Multi-layer perceptron")
legend("topright", legend = layer_names, fill = layer_pal, title = "Layer", cex = 0.8)
```
The degree distribution makes the structural contrast with biological networks immediate:
```{r}
#| label: mlp-vs-celegans-degree
#| fig-cap: "MLP degree distribution (left) vs C. elegans (right). One is a spike, the other a heavy tail. Same type of object — completely different organisation."
p_mlp <- tibble(degree = degree(g_mlp, mode = "total"), network = "MLP (4-6-6-2)")
p_cel <- tibble(degree = degree(giant, mode = "total"), network = "C. elegans")
bind_rows(p_mlp, p_cel) |>
ggplot(aes(x = degree)) +
geom_histogram(bins = 20, fill = "steelblue", colour = "white") +
facet_wrap(~ network, scales = "free") +
labs(x = "Degree", y = "Count",
title = "Degree distributions: artificial vs biological network") +
theme_net()
```
The MLP is **complicated** — many edges, but completely regular. The connectome is **complex** — sparse, heterogeneous, with hubs that matter.
## Message passing: how information flows through a graph
The most important conceptual development linking network science to modern AI is **message passing**. The idea is simple:
1. Each node holds a value (a signal, a feature).
2. At each step, every node collects the values of its neighbours and aggregates them (sum, mean, max).
3. The node updates its own value based on what it received.
4. Repeat for $k$ steps.
After $k$ steps, a node's value reflects information from its $k$-hop neighbourhood. This is not new to biologists — it is how signals propagate through neural circuits, how activation cascades through a signalling network, how a disturbance spreads through a food web.
We can run this directly on the *C. elegans* connectome. Start each neuron with its betweenness centrality. After one step of message passing, each neuron holds the mean betweenness of its direct neighbours.
```{r}
#| label: message-passing
#| fig-cap: "One step of message passing on C. elegans. Each neuron's value before and after receiving neighbour signals. High-betweenness neurons pull their neighbours' values up."
# Initial signal: normalised betweenness
btw <- betweenness(giant, directed = TRUE, normalized = TRUE)
# Adjacency matrix (unweighted, row = sender, col = receiver)
A <- as_adjacency_matrix(giant, sparse = FALSE)
# One message-passing step: each node averages its in-neighbours' signals
deg_in <- degree(giant, mode = "in")
deg_in[deg_in == 0] <- 1 # avoid division by zero for source nodes
msg1 <- as.vector((t(A) %*% btw) / deg_in)
tibble(
neuron = names(btw),
before = as.numeric(btw),
after_1 = msg1
) |>
ggplot(aes(x = before, y = after_1)) +
geom_point(colour = "steelblue", alpha = 0.6) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", colour = "firebrick") +
labs(x = "Betweenness (original)",
y = "Mean neighbour betweenness (after 1 step)",
title = "Message passing on C. elegans — 1 hop") +
theme_net()
```
Neurons above the dashed line are connected to high-betweenness hubs. Neurons below are connected to peripheral nodes. **The graph structure determines what information reaches each node** — independently of any learning.
This is the operation at the heart of **Graph Neural Networks** (GNNs): learn what to aggregate from neighbours, and how to update. Applied to a protein contact graph, this is how AlphaFold 2 predicts protein structure — amino acids are nodes, proximity in 3D space defines edges, message passing over this graph propagates structural context to every residue. The broader programme of applying network thinking to disease and drug targets is network medicine [@Barabasi2015; @Ideker2012a].
## What network science and relational AI share
```{r}
#| label: comparison-table
tribble(
~Concept, ~`Network science`, ~`Graph Neural Networks`,
"Node", "Neuron, protein, species", "Amino acid, atom, gene",
"Edge", "Synapse, interaction, flux", "Contact, bond, co-expression",
"Centrality", "Fixed metric (degree, betw.)", "Learned aggregation function",
"Community", "Dense subgraph partition", "Latent cluster in embedding",
"Information flow", "Shortest paths, diffusion", "Iterated message passing",
"Key question", "What is the structure?", "What does the structure predict?"
) |> kable(caption = "Network science and GNNs: same objects, different questions")
```
The shift is in the question, not the object. Network science asks: *what is the structure of this system?* Relational AI asks: *given this structure, what can we predict?* The graph is the same. The analytical machinery carries over directly.
> **Exercise 5**: Run two steps of message passing instead of one (apply the aggregation twice). Do the values converge or diverge? What does convergence mean biologically — and why might too many steps be a problem for a real GNN on a protein graph?
---
# Epilogue · The Unstable Landscape of Tools {.unnumbered}
In Q1 2026 — the same quarter this workshop was taught — global venture capital investment hit **$297 billion**, of which **$239 billion (81%)** flowed into AI startups. Four of the five largest private funding rounds in history closed in a single 90-day window.
```{r}
#| label: ai-funding
#| fig-cap: "Q1 2026 venture capital: four companies absorbed more than half of all global startup investment. Source: Crunchbase, TechInsider, IndexBox (April 2026)."
#| code-fold: true
tribble(
~Company, ~Billion_USD, ~Category,
"OpenAI", 122, "AI",
"Anthropic", 30, "AI",
"xAI", 20, "AI",
"Waymo", 16, "AI",
"Other AI", 51, "AI",
"Non-AI", 58, "Other"
) |>
mutate(Company = reorder(Company, Billion_USD)) |>
ggplot(aes(x = Company, y = Billion_USD, fill = Category)) +
geom_col() +
geom_text(aes(label = paste0("$", Billion_USD, "B")),
hjust = -0.1, size = 3.5) +
scale_fill_manual(values = c(AI = "steelblue", Other = "grey70")) +
scale_y_continuous(limits = c(0, 145), labels = function(x) paste0("$", x, "B")) +
coord_flip() +
labs(x = NULL, y = "USD (billions)", fill = NULL,
title = "Q1 2026 global VC funding — $297B total",
subtitle = "81% captured by AI; OpenAI alone raised more than all non-AI startups combined") +
theme_net() +
theme(legend.position = "bottom")
```
These numbers are not a celebration. They are a map of concentration and fragility.
**Tools come and go.** The companies absorbing this capital are building platforms, APIs, and proprietary models. They will be acquired, pivoted, deprecated, or superseded — on timescales shorter than a PhD. The history of bioinformatics is a graveyard of tools that were indispensable and then vanished: Galaxy workflows that no longer run, R packages removed from CRAN, web servers that returned 404 before their papers were even cited ten times. AI tools will follow the same trajectory, only faster, because the capital pressure to replace them is now measured in hundreds of billions per quarter.
**What is irreplaceable is data.** Not the model weights, not the API, not the interface — the data. BioGRID, the SNAP connectome library, UniProt, the PDB: open, documented, machine-readable, format-stable datasets that have outlasted dozens of tool generations. The PSI-MI MITAB format we read in Part 4 was defined in 2004 and still works today. The edge list format we used for *C. elegans* is a text file with two columns. It will be readable in 2050.
The key word is *standardised*. A dataset with no metadata is not data — it is noise. The [Genomic Standards Consortium](https://www.gensc.org) (GSC) formalised this insight into the **MIxS** framework: *Minimum Information about any (x) Sequence* [@Yilmaz2011]. Whether the sequence is an amplicon, a metagenome, or a single genome, MIxS defines the minimum metadata fields — sample origin, environment, collection date, sequencing method — that make a record reusable by anyone, with any tool, ten years later. This is technology-agnostic design: the standard does not care whether the sequences were produced on a 454, an Illumina, or a platform that does not exist yet. The metadata schema outlives every instrument generation.
**Reproducibility and open practice** are the institutional expression of the same principle. The [Turing Way](https://the-turing-way.netlify.app) is a community-maintained handbook covering reproducible workflows, open data, version control, and FAIR principles [@TheTuringWay2022]. It is itself an example of what it advocates: free, open-source, collaboratively written, and designed to remain useful regardless of which software stack its readers use.
The practical upshot for your career:
| Invest in | Avoid depending on |
|---|---|
| Standardised, open data formats | Proprietary APIs and cloud platforms |
| Mathematical and biological literacy | Any single tool or framework |
| Version-controlled, documented datasets | Model weights without provenance |
| The concepts behind algorithms | The algorithms' current implementation |
**AI is not a gamble for its users.** It is a gamble taken by a handful of corporations betting hundreds of billions on a single technological trajectory. The risk is theirs. The utility — or the disruption — lands on everyone else. As scientists, the most robust position is not to chase the tools but to produce the kind of data and the kind of rigorous analysis that retains value regardless of which tool reads it next.
The graph formalism we used today — adjacency matrices, edge lists, centralities, communities — predates the computer. It will outlast the current AI cycle. Learn the mathematics. Curate your data. The tools will sort themselves out.
---
# References {.unnumbered}
