---
title: "Composing immune modules with AINet"
author:
  name: Nick Borcherding
  email: ncborch@gmail.com
  affiliation: Washington University in St. Louis, School of Medicine, St. Louis, MO, USA
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc_float: true
package: bHIVE
vignette: >
  %\VignetteIndexEntry{Composing immune modules with AINet}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(message = FALSE, tidy = FALSE)
library(BiocStyle)
```

```{r packages}
library(bHIVE)
```

# Introduction

The bHIVE package has two interfaces. The functional `bHIVE()` runs a whole
analysis in one call, and the introductory vignette,
`vignette("bHIVE", package = "bHIVE")`, uses it on mass cytometry data. The
R6 class `AINet` exposes the same engine as a model object. This vignette
covers `AINet` and the nine immune modules it accepts, each of which
changes one stage of the algorithm. We use the small `iris` data so each
example runs in a second.

Each module is an independent R6 class that encapsulates a single
immunological mechanism. You create modules separately and pass them to
`AINet$new()`, allowing you to mix and match any combination.

```{r module-overview, eval=FALSE}
model <- AINet$new(
  nAntibodies = 20,
  shm            = SHMEngine$new(method = "adaptive"),
  idiotypic      = IdiotypicNetwork$new(),
  germinalCenter = GerminalCenter$new(),
  microenvironment = Microenvironment$new(),
  init           = VDJLibrary$new(method = "pca"),
  activation     = ActivationGate$new(),
  memory         = MemoryPool$new()
)
```

# Somatic Hypermutation Engine

The `SHMEngine` controls how antibodies mutate after cloning. Five
strategies are available, each grounded in a different aspect of SHM biology:

| Method | Biological Analogy | When to Use |
|:-------|:-------------------|:------------|
| `"uniform"` | Random point mutations | Baseline; simple and robust |
| `"airs"` | Affinity-proportional mutation (AIRS) | When high-affinity antibodies should explore less |
| `"hotspot"` | AID targets WRCY motifs | When some features matter more than others |
| `"energy"` | E_SHM ~ N_Mut^2 budget constraint | When you want bounded total perturbation |
| `"adaptive"` | Per-feature Adam-like moment tracking | Best for complex landscapes; most novel |

The **adaptive** method is the most distinctive in which somatic hypermutation and the Adam optimizer function the same. Each feature dimension maintains running mean and variance of past gradients.

```{r shm-comparison, fig.width=8, fig.height=4}
data(iris)
X <- as.matrix(iris[, 1:4])
y <- iris$Species

set.seed(3)
methods <- c("uniform", "airs", "adaptive")
results <- lapply(methods, function(m) {
  shm <- SHMEngine$new(method = m)
  model <- AINet$new(nAntibodies = 15, 
                     maxIter = 15, 
                     shm = shm,
                     verbose = FALSE)
  model$fit(X, y, task = "classification")
  data.frame(
    method   = m,
    accuracy = mean(model$result$assignments == as.character(y)),
    n_ab     = nrow(model$repertoire$as_matrix())
  )
})
do.call(rbind, results)
```

# Idiotypic Network Regulation

In classical AIS, similar antibodies are suppressed via a simple distance
threshold (epsilon). The `IdiotypicNetwork` module replaces this with
network dynamics inspired by Varela and Coutinho's (1991) second-generation
immune network theory.

The key idea is a bell-shaped activation function. Too little stimulation
from other antibodies leads to death, moderate stimulation to activation, and
excessive stimulation to suppression.

The module is experimental. It applies the activation to each pair of
antibodies, so it works as intended on a small, tight repertoire. The example
below builds one: four antibodies at pairwise affinity 0.25, in the middle of
the default activation window, and one far from the rest.

```{r idiotypic-small}
a <- sqrt(log(4) / 2)
A <- rbind(cbind(diag(4) * a, 0), c(0, 0, 0, 0, 10))
rep <- ImmuneRepertoire$new(A)
idi <- IdiotypicNetwork$new(theta_low = 0.01, theta_high = 0.5)
idi$regulate(rep, "gaussian", list(alpha = 1))
idi
round(as.vector(idi$get_population()), 2)
```

The four antibodies that stimulate each other survive. The isolated one is
under-stimulated and dies.

On a full repertoire most pairs of antibodies are far apart. Their affinity
is near zero, below `theta_low`, so each pair counts as suppression and the
dynamics remove every antibody. `AINet` then keeps the top 10% by population
and reports it:

```{r idiotypic-full, message = TRUE}
set.seed(1)
model_idi <- AINet$new(nAntibodies = 25, maxIter = 15,
                       idiotypic = IdiotypicNetwork$new(),
                       verbose = FALSE)
model_idi$fit(X, task = "clustering")
model_idi$repertoire$size()
```

Check that antibodies survive on your data before relying on this module. A
field-based form, in which the activation acts on the total stimulation each
antibody receives, is planned.

# Germinal Center Selection

The `GerminalCenter` module models the competition for T follicular helper
(Tfh) cell help that occurs in real germinal centers. Only antibodies that
receive Tfh help survive -- this implements quality-based selection pressure.

The quality metric is task-aware:

- **Classification**: proportion of correctly matched labels
- **Clustering**: local density (sum of affinities to assigned data points)

```{r gc-demo}
gc <- GerminalCenter$new(
  nTfh = 20,                 # number of Tfh selectors
  selectionPressure = 0.2,   # 0 = no selection, 1 = very strict
  rounds = 1
)

set.seed(4)
model_gc <- AINet$new(
  nAntibodies = 30, 
  maxIter = 15,
  germinalCenter = gc,
  verbose = FALSE
)
model_gc$fit(X, y, task = "classification")
c(antibodies = model_gc$repertoire$size(),
  accuracy = mean(model_gc$predict(X) == y))
```

Selection pressure trades repertoire size for quality. Strong pressure trims
the repertoire hard and, on a small dataset like iris, can cost accuracy.
Tune it against a model without the module.

# Microenvironment

The `Microenvironment` module classifies the local context around each
antibody into three zones and adapts mutation rates accordingly:

- **Stable** (high local density): reduce mutation rate to preserve
  good antibodies
- **Explore** (low local density): increase mutation rate to search
  sparse regions
- **Boundary** (intermediate): moderate mutation, potential trigger for
  class switching

```{r microenv-demo}
me <- Microenvironment$new(
  high_density_threshold = 0.75,
  low_density_threshold = 0.25,
  stabilization_factor = 0.3,   # reduce mutation to 30% in stable zones
  exploration_factor = 2.0      # double mutation in exploration zones
)

rep <- ImmuneRepertoire$new(X[sample(150, 15), ])
env <- me$assess(rep, X)
table(env$zones)
```

# V(D)J Gene Library Initialization

Instead of random sampling, `VDJLibrary` generates antibodies by
combinatorial assembly of gene segments, mimicking V(D)J recombination.
The feature space is split into three segments (V, D, J), each discretized
into alleles, and new antibodies are built by combining one allele from
each segment.

```{r vdj-demo}
vdj <- VDJLibrary$new(nV = 5, 
                      nD = 3, 
                      nJ = 3, 
                      method = "pca")
A <- vdj$generate(20, X)
dim(A)  # 20 antibodies x 4 features
print(vdj)
```

Three decomposition methods are available:

- `"pca"`: PCA components define gene segments (best for correlated features)
- `"cluster"`: k-means within dimension groups creates alleles
- `"random_partition"`: random feature grouping (simplest)

# Two-Signal Activation Gate

In real immunity, B cells require two signals to activate: antigen
recognition (Signal 1) *and* costimulatory context (Signal 2). The
`ActivationGate` module prevents spurious cloning on isolated outliers.

Signal 2 options:

- `"density"`: local data density around the antibody
- `"danger"`: user-provided danger score per data point
- `"entropy"`: local label entropy (classification only)

```{r activation-demo}
gate <- ActivationGate$new(
  signal2_type = "density",
  threshold1 = 0.1,    # minimum affinity (Signal 1)
  threshold2 = 0.3     # density percentile (Signal 2)
)

A <- X[sample(150, 10), ]
rep <- ImmuneRepertoire$new(A)
aff <- rep$affinity_matrix(X, "gaussian")
activated <- gate$evaluate(aff, X, A)
c(activated = sum(activated), total = length(activated))
```

# Memory Pool

The `MemoryPool` archives high-performing antibodies and recalls them
when relevant to current data. This is useful for streaming or
non-stationary data where distribution shifts may occur.

```{r memory-demo}
mp <- MemoryPool$new(
  archive_threshold = 0.01,  # minimum average affinity to archive
  max_memory = 50,
  recall_threshold = 0.01
)

rep <- ImmuneRepertoire$new(X[sample(150, 10), ])
mp$archive(rep, X)
mp

# Recall memories relevant to a subset
recalled <- mp$recall(X[1:30, ])
nrow(recalled)
```

# Class Switcher

The `ClassSwitcher` module changes the effective kernel width of antibodies
based on their microenvironment zone, analogous to real B cell isotype
switching:

| Isotype | Kernel Width | Matching Style |
|:--------|:-------------|:---------------|
| IgM | alpha = 0.1 | Broad (exploration) |
| IgG | alpha = 5.0 | Specific (discrimination) |
| IgA | alpha = 1.0 | Intermediate (boundary patrol) |

```{r classswitcher-demo}
cs <- ClassSwitcher$new(alpha_IgM = 0.1, alpha_IgG = 5.0, alpha_IgA = 1.0)

rep <- ImmuneRepertoire$new(X[sample(150, 10), ])
zones <- sample(c("stable", "explore", "boundary"), 10, replace = TRUE)
alphas <- cs$switch_isotypes(rep, zones)
data.frame(zone = zones, isotype = rep$metadata$isotype, alpha = alphas)
```

# Convergent Selection

The `ConvergentSelector` identifies "public antibodies" -- antibodies that
appear across multiple independent bHIVE runs. This implements the
immunological concept of public clonotypes as a biologically-motivated
ensemble method.

```{r convergent-demo}
# Run 3 independent bHIVE analyses
set.seed(42)
results <- lapply(1:3, function(i) {
  m <- AINet$new(nAntibodies = 15, maxIter = 10, verbose = FALSE)
  m$fit(X, task = "clustering")
  m$result
})

conv <- ConvergentSelector$new(tolerance = 1.0, min_appearances = 2)
public <- conv$from_results(results)
nrow(public)
```

# Putting It All Together

Here is a full example combining multiple modules for a classification task:

```{r full-example}
set.seed(2)
model <- AINet$new(
  nAntibodies = 25,
  maxIter = 20,
  k = 3,
  beta = 5,
  shm = SHMEngine$new(method = "adaptive", base_rate = 0.1),
  germinalCenter = GerminalCenter$new(nTfh = 8, selectionPressure = 0.4),
  microenvironment = Microenvironment$new(),
  classSwitcher = ClassSwitcher$new(),
  verbose = FALSE
)

model$fit(X, y, task = "classification")
model
mean(model$predict(X) == y)
```

Not every module combination improves every dataset. Compare a combination
against the plain `AINet` before keeping it, on held-out data or with
cross-validation. `caret::train(method = bHIVEmodel)` passes extra arguments,
modules included, on to `bHIVE()`.

# Session Information

```{r session}
sessionInfo()
```