Hands-On Labs — Building DOGMA from Scratch

Theory without practice is inert. This chapter converts the preceding seventeen chapters of conceptual framework into running code through eight self-contained laboratory exercises. Each lab targets a specific layer of the DOGMA architecture, progressing from the molecular foundations of Watson-Crick base pairing through strand displacement kinetics, chemical reaction networks, tetrahedral memory, thermodynamic attention, and culminating in a complete end-to-end DOGMA system trained on a genomic task. By the final lab, you will have built, trained, visualized, and benchmarked a working DOGMA model from scratch.

Every lab follows a consistent format: learning objectives, step-by-step instructions with complete code, expected outputs that let you verify correctness, and discussion questions that connect implementation to theory. All code is written in Python with PyTorch and runs on CPU.

18.1 Environment Setup

All labs require Python 3.11 or later. Create an isolated virtual environment and install the core dependencies:

1# Environment setup
2python3.11 -m venv dogma-env && source dogma-env/bin/activate
3pip install numpy>=1.25 scipy>=1.11 torch>=2.1
4pip install matplotlib seaborn scikit-learn
5pip install dataclasses-json pyyaml pytest
6
7# Recommended project structure
8# dogma-labs/
9#   lab1_base_pairing.py
10#   lab2_strand_displacement.py
11#   lab3_crn_logic_gate.py
12#   lab4_tetra_memory.py
13#   lab5_thermo_attention.py
14#   lab6_dogma_pipeline.py
15#   lab7_dogma_vs_transformer.py
16#   lab8_visualization.py
17#   utils.py
18#   data/           <- toy datasets for Labs 6-8
19#   figures/        <- saved plots

18.2 Lab 1: Implementing Watson-Crick Base Pairing in Code

18.2.1 Learning Objectives

Upon completing this lab, you will be able to:

1.

Encode arbitrary DNA sequences as numerical tensors using one-hot and embedding representations.

2.

Compute the Watson-Crick complement of any DNA strand programmatically.

3.

Simulate hybridization between two strands using a dot-product affinity model.

4.

Quantify hybridization strength and identify mismatches.

18.2.2 Background

Watson-Crick base pairing is the foundation of all DNA computation (Chapter 1). Adenine (A) pairs with Thymine (T), and Cytosine (C) pairs with Guanine (G). Each pairing is mediated by hydrogen bonds—two for A-T and three for C-G—giving C-G pairs higher thermodynamic stability. In DOGMA, this complementarity principle is abstracted into the affinity function that determines which genomic bases interact.

18.2.3 Step-by-Step Instructions

Step 1: Define the encoding scheme. We represent each nucleotide as a one-hot vector in \(\mathbb {R}^4\) and define the complement mapping.

1import torch
2import torch.nn.functional as F
3import numpy as np
4
5# Nucleotide-to-index mapping
6NUC_TO_IDX = {’A’: 0, ’T’: 1, ’C’: 2, ’G’: 3}
7IDX_TO_NUC = {v: k for k, v in NUC_TO_IDX.items()}
8
9# Watson-Crick complement mapping
10WC_COMPLEMENT = {’A’: ’T’, ’T’: ’A’, ’C’: ’G’, ’G’: ’C’}
11
12def encode_sequence(seq: str) -> torch.Tensor:
13    """Encode a DNA string as one-hot vectors.
14
15    Args:
16        seq: DNA string, e.g. "ATCGATCG"
17
18    Returns:
19        Tensor of shape (L, 4) where L = len(seq)
20    """
21    indices = [NUC_TO_IDX[c] for c in seq.upper()]
22    return F.one_hot(torch.tensor(indices), num_classes=4).float()
23
24def decode_sequence(encoded: torch.Tensor) -> str:
25    """Convert one-hot tensor back to DNA string."""
26    indices = encoded.argmax(dim=-1).tolist()
27    return ’’.join(IDX_TO_NUC[i] for i in indices)
28
29def complement(seq: str) -> str:
30    """Return the Watson-Crick complement of a DNA string."""
31    return ’’.join(WC_COMPLEMENT[c] for c in seq.upper())
32
33def reverse_complement(seq: str) -> str:
34    """Return the reverse complement (antiparallel strand)."""
35    return complement(seq)[::-1]

Step 2: Build the complementarity matrix. The Watson-Crick pairing rules can be expressed as a \(4 \times 4\) matrix \(\mathbf {W}\) where \(W_{ij} = 1\) if nucleotides \(i\) and \(j\) are complementary and \(0\) otherwise:

\begin {equation} \mathbf {W} = \begin {pmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end {pmatrix}, \quad \text {row/col order: A, T, C, G}. \label {eq:wc-matrix} \end {equation}

1# Complementarity matrix (Eq. 18.1)
2WC_MATRIX = torch.tensor([
3    [0, 1, 0, 0],  # A pairs with T
4    [1, 0, 0, 0],  # T pairs with A
5    [0, 0, 0, 1],  # C pairs with G
6    [0, 0, 1, 0],  # G pairs with C
7], dtype=torch.float32)
8
9def hybridization_affinity(seq1: str, seq2: str) -> torch.Tensor:
10    """Compute per-position hybridization affinity.
11
12    For each position i, affinity[i] = 1.0 if seq1[i] and
13    seq2[i] are Watson-Crick complements, else 0.0.
14
15    Args:
16        seq1: First DNA strand  (5’ -> 3’)
17        seq2: Second DNA strand (3’ -> 5’, i.e., antiparallel)
18
19    Returns:
20        Tensor of shape (L,) with per-position affinities
21    """
22    enc1 = encode_sequence(seq1)   # (L, 4)
23    enc2 = encode_sequence(seq2)   # (L, 4)
24    # For each position: enc1[i] @ W @ enc2[i]^T
25    affinity = torch.einsum(’ij,jk,ik->i’, enc1, WC_MATRIX, enc2)
26    return affinity
27
28def hybridization_score(seq1: str, seq2: str) -> float:
29    """Overall hybridization score = fraction of matched pairs."""
30    aff = hybridization_affinity(seq1, seq2)
31    return aff.mean().item()
32
33def find_mismatches(seq1: str, seq2: str) -> list[int]:
34    """Return indices where base pairing fails."""
35    aff = hybridization_affinity(seq1, seq2)
36    return (aff == 0).nonzero(as_tuple=True)[0].tolist()

Step 3: Test hybridization.

1# Perfect complement
2strand_a = "ATCGATCG"
3strand_b = reverse_complement(strand_a)
4print(f"Strand A:    5’-{strand_a}-3’")
5print(f"Strand B:    3’-{strand_b}-5’")
6
7score = hybridization_score(strand_a, strand_b[::-1])
8print(f"Hybridization score: {score:.2f}")  # Expected: 1.00
9
10# Introduce a mismatch at position 3
11strand_c = list(strand_b[::-1])
12strand_c[3] = ’A’  # was ’C’, now mismatch
13strand_c = ’’.join(strand_c)
14score_mm = hybridization_score(strand_a, strand_c)
15mismatches = find_mismatches(strand_a, strand_c)
16print(f"With mismatch: score={score_mm:.2f}, "
17      f"mismatch positions={mismatches}")
18# Expected: score=0.88, mismatch positions=[3]

18.2.4 Expected Output

Strand A:    5’-ATCGATCG-3’
Strand B:    3’-TAGCTAGC-5’
Hybridization score: 1.00
With mismatch: score=0.88, mismatch positions=[3]

18.2.5 Discussion Questions

1.

Why does DOGMA use a complementarity matrix rather than a simple lookup table? What advantage does the matrix formulation offer for batch computation?

2.

In real DNA, C-G pairs are stronger than A-T pairs due to three vs. two hydrogen bonds. Modify the complementarity matrix \(\mathbf {W}\) so that C-G matches have weight 1.5 and A-T matches have weight 1.0. How does this change the hybridization score for GC-rich vs. AT-rich sequences?

3.

How does the reverse_complement function relate to the antiparallel orientation of real DNA strands?

18.3 Lab 2: Building a Toehold-Mediated Strand Displacement Simulator

18.3.1 Learning Objectives

Upon completing this lab, you will be able to:

1.

Explain the three phases of toehold-mediated strand displacement: toehold binding, branch migration, and completion.

2.

Implement a kinetic simulator for strand displacement using ordinary differential equations.

3.

Plot concentration trajectories showing how an invader strand displaces an incumbent.

4.

Analyze how toehold length affects displacement rate.

18.3.2 Background

Toehold-mediated strand displacement (Chapter 2) is the mechanism by which one DNA strand replaces another on a complementary template. The process proceeds in three phases:

1.

Toehold binding. The invader strand binds to a short single-stranded region (the toehold, typically 5–8 nucleotides) on the target complex. This initiates the reaction.

2.

Branch migration. The invader progressively displaces the incumbent strand through a random walk of base-pair exchanges. Each step is reversible, but the toehold provides the thermodynamic bias that drives the reaction forward.

3.

Completion. The incumbent strand is fully displaced and released into solution. The invader is now hybridized to the template.

We model this process using mass-action kinetics. Let \([S]\) denote the substrate (template:incumbent complex), \([I]\) the invader concentration, \([SI]\) the intermediate (toehold-bound state), and \([P]\) the product (template:invader complex plus free incumbent):

\begin {align} S + I &\xrightarrow {k_{\text {bind}}} SI, \label {eq:sd-bind} \\ SI &\xrightarrow {k_{\text {migrate}}} P. \label {eq:sd-migrate} \end {align}

18.3.3 Step-by-Step Instructions

Step 1: Define the ODE system.

1import numpy as np
2from scipy.integrate import solve_ivp
3import matplotlib.pyplot as plt
4
5def strand_displacement_odes(t, y, k_bind, k_migrate):
6    """ODE system for toehold-mediated strand displacement.
7
8    Species: y = [S, I, SI, P]
9      S  = substrate (template:incumbent)
10      I  = invader strand
11      SI = intermediate (toehold-bound)
12      P  = product (template:invader + free incumbent)
13    """
14    S, I, SI, P = y
15    # Toehold binding: S + I -> SI
16    r_bind = k_bind * S * I
17    # Branch migration: SI -> P
18    r_migrate = k_migrate * SI
19
20    dS  = -r_bind
21    dI  = -r_bind
22    dSI = r_bind - r_migrate
23    dP  = r_migrate
24    return [dS, dI, dSI, dP]

Step 2: Simulate and plot the kinetics.

1def simulate_displacement(k_bind, k_migrate, S0=1.0, I0=1.5,
2                          t_end=100.0, n_points=500):
3    """Simulate strand displacement and return trajectories."""
4    y0 = [S0, I0, 0.0, 0.0]  # initial concentrations
5    t_span = (0, t_end)
6    t_eval = np.linspace(0, t_end, n_points)
7
8    sol = solve_ivp(
9        strand_displacement_odes, t_span, y0,
10        args=(k_bind, k_migrate),
11        t_eval=t_eval, method=’RK45’
12    )
13    return sol
14
15def plot_displacement(sol, title="Strand Displacement Kinetics"):
16    """Plot concentration trajectories."""
17    labels = [’Substrate [S]’, ’Invader [I]’,
18              ’Intermediate [SI]’, ’Product [P]’]
19    colors = [’#2196F3’, ’#FF9800’, ’#4CAF50’, ’#F44336’]
20
21    fig, ax = plt.subplots(figsize=(8, 5))
22    for i, (label, color) in enumerate(zip(labels, colors)):
23        ax.plot(sol.t, sol.y[i], label=label,
24                color=color, linewidth=2)
25    ax.set_xlabel(’Time (arbitrary units)’, fontsize=12)
26    ax.set_ylabel(’Concentration (nM)’, fontsize=12)
27    ax.set_title(title, fontsize=14)
28    ax.legend(fontsize=11)
29    ax.grid(True, alpha=0.3)
30    plt.tight_layout()
31    plt.savefig(’figures/lab2_displacement.png’, dpi=150)
32    plt.show()
33
34# Run simulation with default parameters
35sol = simulate_displacement(k_bind=0.05, k_migrate=0.02)
36plot_displacement(sol)

Step 3: Analyze toehold length dependence. The binding rate \(k_{\text {bind}}\) increases approximately exponentially with toehold length \(n\), following \(k_{\text {bind}} \approx k_0 \cdot 3^{n}\) for \(n \in [1, 7]\) (Chapter 2):

1def toehold_sweep(toehold_lengths=range(1, 9),
2                  k0=1e-4, k_migrate=0.02):
3    """Sweep toehold length, measure time to 90% completion."""
4    half_times = []
5    for n in toehold_lengths:
6        k_bind = k0 * (3 ** n)
7        sol = simulate_displacement(k_bind, k_migrate, t_end=500)
8        # Find time to reach 90% product
9        target = 0.9 * min(sol.y[0][0], sol.y[1][0])
10        idx = np.searchsorted(sol.y[3], target)
11        t90 = sol.t[min(idx, len(sol.t) - 1)]
12        half_times.append(t90)
13        print(f"Toehold {n}nt: k_bind={k_bind:.4f}, "
14              f"t_90={t90:.1f}")
15
16    fig, ax = plt.subplots(figsize=(7, 4))
17    ax.plot(list(toehold_lengths), half_times, ’o-’,
18            color=’#2196F3’, linewidth=2, markersize=8)
19    ax.set_xlabel(’Toehold Length (nt)’, fontsize=12)
20    ax.set_ylabel(’Time to 90%% Completion’, fontsize=12)
21    ax.set_title(’Displacement Rate vs Toehold Length’, fontsize=14)
22    ax.set_yscale(’log’)
23    ax.grid(True, alpha=0.3)
24    plt.tight_layout()
25    plt.savefig(’figures/lab2_toehold_sweep.png’, dpi=150)
26    plt.show()
27
28toehold_sweep()

18.3.4 Expected Output

The concentration plot should show: (1) Substrate and Invader concentrations decreasing together as they bind; (2) Intermediate concentration rising then falling as the toehold-bound complex forms and then undergoes branch migration; (3) Product concentration rising in a sigmoidal curve to its final value. The toehold sweep should show that the time to 90% completion drops exponentially as toehold length increases from 1 to 7 nucleotides, then plateaus beyond 7 nt.

18.3.5 Discussion Questions

1.

Why is the intermediate \([SI]\) transient—why does it rise and then fall? Relate this to the concept of a kinetic bottleneck.

2.

In DOGMA’s regulation engine, how does toehold-mediated displacement inspire the mechanism by which a new context signal “displaces” an old regulatory state?

3.

Modify the model to include a reverse reaction \(SI \xrightarrow {k_{\text {rev}}} S + I\). How does the equilibrium change?

4.

Real branch migration is a random walk. Why is a deterministic ODE model still useful?

18.4 Lab 3: Implementing a DNA Logic Gate (AND) as a CRN

18.4.1 Learning Objectives

Upon completing this lab, you will be able to:

1.

Translate a Boolean AND gate into a set of chemical reactions with mass-action kinetics.

2.

Simulate the chemical reaction network (CRN) using numerical ODE integration.

3.

Verify that the CRN correctly computes the AND truth table for all four input combinations.

4.

Explain how CRN-based logic connects to DOGMA’s regulatory logic.

18.4.2 Background

Chemical reaction networks (CRNs) provide a Turing-complete model of computation (Chapter 2). A DNA AND gate can be realized with two input species \(X_1\) and \(X_2\) and one output species \(Y\), mediated by a catalyst \(C\):

\begin {align} X_1 + X_2 &\xrightarrow {k_1} C, \label {eq:and-step1} \\ C &\xrightarrow {k_2} Y + C. \label {eq:and-step2} \end {align}

The first reaction requires both inputs to be present (implementing the AND logic). The second reaction catalytically produces the output. When either input is absent (concentration \(\approx 0\)), the first reaction cannot proceed, and \(Y\) remains at zero.

18.4.3 Step-by-Step Instructions

Step 1: Define the CRN ODE system.

1import numpy as np
2from scipy.integrate import solve_ivp
3import matplotlib.pyplot as plt
4
5def and_gate_odes(t, y, k1, k2):
6    """CRN for a DNA AND gate.
7
8    Species: y = [X1, X2, C, Y]
9      X1, X2 = input signals
10      C      = intermediate catalyst
11      Y      = output signal
12    """
13    X1, X2, C, Y = y
14    # Reaction 1: X1 + X2 -> C
15    r1 = k1 * X1 * X2
16    # Reaction 2: C -> Y + C  (catalytic)
17    r2 = k2 * C
18
19    dX1 = -r1
20    dX2 = -r1
21    dC  = r1          # C is produced in r1, not consumed in r2
22    dY  = r2
23    return [dX1, dX2, dC, dY]
24
25def simulate_and_gate(X1_init, X2_init, k1=0.1, k2=0.05,
26                      t_end=200, n_points=500):
27    """Simulate the AND gate for given input concentrations."""
28    y0 = [X1_init, X2_init, 0.0, 0.0]
29    t_eval = np.linspace(0, t_end, n_points)
30    sol = solve_ivp(and_gate_odes, (0, t_end), y0,
31                    args=(k1, k2), t_eval=t_eval, method=’RK45’)
32    return sol

Step 2: Verify the truth table.

1def verify_truth_table(threshold=0.3):
2    """Verify the AND gate for all input combinations."""
3    inputs = [(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0)]
4    expected = [0, 0, 0, 1]
5
6    fig, axes = plt.subplots(2, 2, figsize=(12, 9))
7    axes = axes.flatten()
8
9    print("AND Gate Truth Table Verification")
10    print("-" * 45)
11    print(f"{’X1’:>4} {’X2’:>4} {’Y_final’:>10} "
12          f"{’Expected’:>10} {’Pass’:>6}")
13    print("-" * 45)
14
15    for idx, ((x1, x2), exp) in enumerate(
16            zip(inputs, expected)):
17        sol = simulate_and_gate(x1, x2)
18        y_final = sol.y[3][-1]
19        digital = 1 if y_final > threshold else 0
20        passed = digital == exp
21
22        print(f"{x1:4.1f} {x2:4.1f} {y_final:10.4f} "
23              f"{exp:10d} {’OK’ if passed else ’FAIL’:>6}")
24
25        ax = axes[idx]
26        labels = [’$X_1$’, ’$X_2$’, ’$C$’, ’$Y$’]
27        colors = [’#2196F3’, ’#FF9800’, ’#9C27B0’, ’#F44336’]
28        for i, (lbl, clr) in enumerate(zip(labels, colors)):
29            ax.plot(sol.t, sol.y[i], label=lbl,
30                    color=clr, linewidth=2)
31        ax.set_title(f’$X_1$={x1:.0f}, $X_2$={x2:.0f}’,
32                     fontsize=12)
33        ax.set_xlabel(’Time’)
34        ax.set_ylabel(’Concentration’)
35        ax.legend(fontsize=9)
36        ax.grid(True, alpha=0.3)
37
38    plt.suptitle(’AND Gate CRN: All Input Combinations’,
39                 fontsize=14, y=1.02)
40    plt.tight_layout()
41    plt.savefig(’figures/lab3_and_gate.png’, dpi=150)
42    plt.show()
43
44verify_truth_table()

18.4.4 Expected Output

AND Gate Truth Table Verification
---------------------------------------------
  X1   X2    Y_final   Expected   Pass
---------------------------------------------
 0.0  0.0     0.0000          0     OK
 1.0  0.0     0.0000          0     OK
 0.0  1.0     0.0000          0     OK
 1.0  1.0     0.8347          1     OK

The four subplots should show: (1) all species flat at zero for (0,0); (2–3) only one input present, no output produced; (4) both inputs consumed, catalyst rising, output \(Y\) growing.

18.4.5 Discussion Questions

1.

Modify the CRN to implement an OR gate. Which reaction needs to change?

2.

The second reaction in the AND-gate system amplifies the signal. What happens if you remove the catalytic cycle and make \(C \to Y\) a stoichiometric conversion instead?

3.

In DOGMA, the regulation engine applies multiple regulatory elements to compute an activation map. How is this analogous to a CRN where multiple “input species” determine whether a “gene” (output species) is activated?

4.

Real DNA logic gates suffer from signal leakage—small amounts of output even when inputs are nominally zero. Add a leakage term \(r_{\text {leak}} = k_{\text {leak}} \cdot 0.01\) to the ODE for \(Y\) and observe how it affects the truth table.

18.5 Lab 4: Building TetraMemory from Scratch

18.5.1 Learning Objectives

Upon completing this lab, you will be able to:

1.

Construct the \(K_4\)-complete graph data structure that underlies TetraMemory.

2.

Implement read and write operations on tetrahedral memory cells.

3.

Demonstrate that six edges in a 4-vertex complete graph store six independent values.

4.

Perform associative retrieval by querying a single vertex to obtain its three connected edges.

18.5.2 Background

TetraMemory (Chapter 10) uses the complete graph \(K_4\) as its fundamental memory unit. A tetrahedron has four vertices and six edges. Each vertex stores a key embedding, and each edge stores a relational value. This gives \(\binom {4}{2} = 6\) independent storage slots per tetrahedron, making it more efficient than flat key-value storage for relational data.

The key insight is that reading from one vertex retrieves three edges simultaneously—the associative fan-out property. Writing to an edge updates the relationship between exactly two vertices without disturbing the other four edges.

18.5.3 Step-by-Step Instructions

Step 1: Implement the \(K_4\) data structure.

1import torch
2import torch.nn.functional as F
3from dataclasses import dataclass, field
4from itertools import combinations
5
6@dataclass
7class TetraCell:
8    """A single tetrahedral memory cell (K4 complete graph).
9
10    Attributes:
11        vertices: dict mapping vertex ID (0-3) to key
12                  embedding of shape (d,)
13        edges:    dict mapping edge tuple (i,j) to value
14                  embedding of shape (d,)
15        dim:      embedding dimension
16    """
17    dim: int = 64
18    vertices: dict = field(default_factory=dict)
19    edges: dict = field(default_factory=dict)
20
21    def __post_init__(self):
22        # Initialize 4 vertices with zero embeddings
23        for v in range(4):
24            self.vertices[v] = torch.zeros(self.dim)
25        # Initialize 6 edges (all pairs) with zero embeddings
26        for i, j in combinations(range(4), 2):
27            self.edges[(i, j)] = torch.zeros(self.dim)
28
29    def write_vertex(self, vertex_id: int,
30                     key: torch.Tensor) -> None:
31        """Write a key embedding to a vertex."""
32        assert 0 <= vertex_id <= 3, "Vertex ID must be 0-3"
33        assert key.shape == (self.dim,)
34        self.vertices[vertex_id] = key.clone()
35
36    def write_edge(self, v1: int, v2: int,
37                   value: torch.Tensor) -> None:
38        """Write a value embedding to the edge (v1, v2)."""
39        edge_key = (min(v1, v2), max(v1, v2))
40        assert edge_key in self.edges, f"Invalid edge {edge_key}"
41        assert value.shape == (self.dim,)
42        self.edges[edge_key] = value.clone()
43
44    def read_vertex(self, vertex_id: int) -> dict:
45        """Read a vertex and its three connected edges.
46
47        Returns:
48            dict with ’key’ and ’neighbors’: list of
49            (neighbor_id, edge_value) tuples
50        """
51        key = self.vertices[vertex_id]
52        neighbors = []
53        for other in range(4):
54            if other == vertex_id:
55                continue
56            edge_key = (min(vertex_id, other),
57                        max(vertex_id, other))
58            neighbors.append((other, self.edges[edge_key]))
59        return {’key’: key, ’neighbors’: neighbors}
60
61    def read_edge(self, v1: int, v2: int) -> torch.Tensor:
62        """Read the value stored on edge (v1, v2)."""
63        edge_key = (min(v1, v2), max(v1, v2))
64        return self.edges[edge_key]

Step 2: Implement associative retrieval.

1@dataclass
2class TetraMemory:
3    """A collection of TetraCells with content-addressed access."""
4    dim: int = 64
5    cells: list = field(default_factory=list)
6
7    def add_cell(self) -> TetraCell:
8        """Create and register a new tetrahedral cell."""
9        cell = TetraCell(dim=self.dim)
10        self.cells.append(cell)
11        return cell
12
13    def query(self, query_vec: torch.Tensor,
14              top_k: int = 1) -> list[dict]:
15        """Content-addressed query: find cells whose vertices
16        best match the query vector.
17
18        Returns list of (cell_idx, vertex_id, similarity, data).
19        """
20        results = []
21        for c_idx, cell in enumerate(self.cells):
22            for v_id, key in cell.vertices.items():
23                if key.norm() < 1e-8:
24                    continue  # skip unwritten vertices
25                sim = F.cosine_similarity(
26                    query_vec.unsqueeze(0),
27                    key.unsqueeze(0)
28                ).item()
29                data = cell.read_vertex(v_id)
30                results.append({
31                    ’cell’: c_idx, ’vertex’: v_id,
32                    ’similarity’: sim, ’data’: data
33                })
34        results.sort(key=lambda x: x[’similarity’],
35                     reverse=True)
36        return results[:top_k]
37
38    def utilization(self) -> dict:
39        """Report memory utilization statistics."""
40        total_vertices = len(self.cells) * 4
41        total_edges = len(self.cells) * 6
42        used_v = sum(
43            1 for c in self.cells
44            for v in c.vertices.values()
45            if v.norm() > 1e-8
46        )
47        used_e = sum(
48            1 for c in self.cells
49            for e in c.edges.values()
50            if e.norm() > 1e-8
51        )
52        return {
53            ’cells’: len(self.cells),
54            ’vertices_used’: f"{used_v}/{total_vertices}",
55            ’edges_used’: f"{used_e}/{total_edges}",
56            ’vertex_utilization’: used_v / max(total_vertices, 1),
57            ’edge_utilization’: used_e / max(total_edges, 1),
58        }

Step 3: Test read/write operations.

1torch.manual_seed(42)
2
3# Create memory and a cell
4memory = TetraMemory(dim=64)
5cell = memory.add_cell()
6
7# Write key embeddings to vertices (representing concepts)
8concepts = {
9    0: "protein",   1: "DNA",
10    2: "ribosome",  3: "mRNA"
11}
12embeddings = {}
13for v_id, name in concepts.items():
14    # Use deterministic pseudo-embeddings
15    emb = torch.randn(64)
16    emb = emb / emb.norm()  # normalize
17    cell.write_vertex(v_id, emb)
18    embeddings[name] = emb
19    print(f"Wrote vertex {v_id}: {name}")
20
21# Write relational values to edges
22relations = {
23    (0, 1): "encodes",      (0, 2): "translated_by",
24    (0, 3): "transcribed_from",
25    (1, 2): "read_by",      (1, 3): "transcribed_to",
26    (2, 3): "translates",
27}
28for (v1, v2), rel_name in relations.items():
29    rel_emb = torch.randn(64)
30    rel_emb = rel_emb / rel_emb.norm()
31    cell.write_edge(v1, v2, rel_emb)
32    print(f"Wrote edge ({v1},{v2}): {rel_name}")
33
34# Associative read: query vertex 0 (protein)
35result = cell.read_vertex(0)
36print(f"\nReading vertex 0 (protein):")
37print(f"  Key norm: {result[’key’].norm():.4f}")
38print(f"  Connected to {len(result[’neighbors’])} neighbors")
39for nbr_id, edge_val in result[’neighbors’]:
40    print(f"  -> vertex {nbr_id} ({concepts[nbr_id]}), "
41          f"edge norm: {edge_val.norm():.4f}")
42
43# Content-addressed query
44query = embeddings["DNA"] + 0.1 * torch.randn(64)
45matches = memory.query(query, top_k=2)
46print(f"\nQuery similar to ’DNA’:")
47for m in matches:
48    print(f"  Cell {m[’cell’]}, vertex {m[’vertex’]} "
49          f"({concepts[m[’vertex’]]}), "
50          f"sim={m[’similarity’]:.4f}")
51
52print(f"\nUtilization: {memory.utilization()}")

18.5.4 Expected Output

Wrote vertex 0: protein
Wrote vertex 1: DNA
Wrote vertex 2: ribosome
Wrote vertex 3: mRNA
Wrote edge (0,1): encodes
...
Reading vertex 0 (protein):
  Key norm: 1.0000
  Connected to 3 neighbors
  -> vertex 1 (DNA), edge norm: 1.0000
  -> vertex 2 (ribosome), edge norm: 1.0000
  -> vertex 3 (mRNA), edge norm: 1.0000

Query similar to ’DNA’:
  Cell 0, vertex 1 (DNA), sim=0.9876
  Cell 0, vertex 3 (mRNA), sim=0.1234

Utilization: {’cells’: 1, ’vertices_used’: ’4/4’,
  ’edges_used’: ’6/6’, ’vertex_utilization’: 1.0,
  ’edge_utilization’: 1.0}

18.5.5 Discussion Questions

1.

A single tetrahedron stores 6 edge values. Compare this to a flat key-value store with 4 keys. What relational information is lost in the flat representation?

2.

How does the associative fan-out (reading one vertex retrieves three edges) relate to the “spreading activation” models used in cognitive science?

3.

Extend the TetraCell class to support weighted vertices, where each vertex has a “relevance” score \(r \in [0, 1]\). How would you use this during retrieval?

4.

If you have \(N\) concepts to store, how many tetrahedra do you need, and how should concepts be assigned to vertices to maximize edge utility?

18.6 Lab 5: Implementing the Thermodynamic Attention Mechanism

18.6.1 Learning Objectives

Upon completing this lab, you will be able to:

1.

Implement the SantaLucia nearest-neighbor model for DNA duplex thermodynamics.

2.

Compute \(\Delta G\), \(\Delta H\), and \(\Delta S\) for arbitrary DNA sequences.

3.

Build a thermodynamic attention layer that uses free energy as the attention score.

4.

Compare thermodynamic attention to standard scaled dot-product attention.

18.6.2 Background

Standard transformer attention computes scores as \(\text {softmax}(QK^\top / \sqrt {d_k})\). DOGMA replaces this with a thermodynamically motivated score based on DNA hybridization free energy (Chapter 9). The SantaLucia nearest-neighbor model predicts the Gibbs free energy of duplex formation from dinucleotide parameters: \begin {equation} \Delta G^\circ _{37} = \sum _{i=1}^{n-1} \Delta G^\circ (s_i s_{i+1} / s'_i s'_{i+1}) + \Delta G^\circ _{\text {init}}, \label {eq:ch18-santalucia} \end {equation} where \(s_i s_{i+1}\) is a dinucleotide and \(s'_i s'_{i+1}\) is its complement. The initiation term \(\Delta G^\circ _{\text {init}}\) accounts for the entropic cost of bringing two strands together. More negative \(\Delta G\) means stronger binding and thus higher attention.

18.6.3 Step-by-Step Instructions

Step 1: Encode the SantaLucia nearest-neighbor parameters.

1import torch
2import torch.nn as nn
3import torch.nn.functional as F
4
5# SantaLucia (1998) unified nearest-neighbor parameters
6# Format: ’XY/X’Y’’ -> (dH in kcal/mol, dS in cal/mol/K)
7# where XY is the top strand 5’->3’ dinucleotide
8# and X’Y’ is the bottom strand 3’->5’ complement
9NN_PARAMS = {
10    ’AA/TT’: (-7.9, -22.2),  ’AT/TA’: (-7.2, -20.4),
11    ’TA/AT’: (-7.2, -21.3),  ’CA/GT’: (-8.5, -22.7),
12    ’GT/CA’: (-8.4, -22.4),  ’CT/GA’: (-7.8, -21.0),
13    ’GA/CT’: (-8.2, -22.2),  ’CG/GC’: (-10.6, -27.2),
14    ’GC/CG’: (-9.8, -24.4),  ’GG/CC’: (-8.0, -19.9),
15    ’AC/TG’: (-7.8, -21.0),  ’TC/AG’: (-8.2, -22.2),
16    ’AG/TC’: (-7.8, -21.0),  ’TG/AC’: (-8.5, -22.7),
17    ’TT/AA’: (-7.9, -22.2),  ’CC/GG’: (-8.0, -19.9),
18}
19
20# Initiation parameters
21DG_INIT = 1.96  # kcal/mol (initiation penalty)
22
23WC = {’A’: ’T’, ’T’: ’A’, ’C’: ’G’, ’G’: ’C’}
24
25def compute_thermodynamics(seq: str,
26                           temp_celsius: float = 37.0):
27    """Compute dG, dH, dS for a DNA duplex using the
28    SantaLucia nearest-neighbor model.
29
30    Args:
31        seq: DNA sequence (top strand, 5’->3’)
32        temp_celsius: Temperature in Celsius
33
34    Returns:
35        dict with ’dG’, ’dH’, ’dS’ values
36    """
37    seq = seq.upper()
38    comp = ’’.join(WC[c] for c in seq)
39    temp_k = temp_celsius + 273.15
40
41    total_dh = 0.0  # kcal/mol
42    total_ds = 0.0  # cal/mol/K
43
44    for i in range(len(seq) - 1):
45        dinuc = seq[i:i+2]
46        comp_dinuc = comp[i:i+2]
47        key = f"{dinuc}/{comp_dinuc}"
48        if key in NN_PARAMS:
49            dh, ds = NN_PARAMS[key]
50        else:
51            # Try reverse complement lookup
52            rev_key = f"{comp_dinuc[::-1]}/{dinuc[::-1]}"
53            dh, ds = NN_PARAMS.get(rev_key, (-7.5, -21.0))
54        total_dh += dh
55        total_ds += ds
56
57    # Add initiation
58    total_dg = total_dh - temp_k * (total_ds / 1000.0)
59    total_dg += DG_INIT
60
61    return {
62        ’dG’: total_dg,   # kcal/mol
63        ’dH’: total_dh,   # kcal/mol
64        ’dS’: total_ds,   # cal/mol/K
65        ’Tm’: (total_dh * 1000.0) / total_ds - 273.15
66              if total_ds != 0 else 0.0
67    }

Step 2: Build the thermodynamic attention layer.

1class ThermodynamicAttention(nn.Module):
2    """Attention mechanism using free-energy-inspired scores.
3
4    Instead of dot-product attention, we compute scores
5    based on a learned "hybridization energy" that mimics
6    the SantaLucia nearest-neighbor model.
7
8    Args:
9        d_model: embedding dimension
10        n_heads: number of attention heads
11        temp:    softmax temperature (analogous to
12                 physical temperature)
13    """
14    def __init__(self, d_model: int = 64, n_heads: int = 4,
15                 temp: float = 1.0):
16        super().__init__()
17        assert d_model % n_heads == 0
18        self.d_model = d_model
19        self.n_heads = n_heads
20        self.d_k = d_model // n_heads
21        self.temp = temp
22
23        # Q, K, V projections
24        self.W_q = nn.Linear(d_model, d_model, bias=False)
25        self.W_k = nn.Linear(d_model, d_model, bias=False)
26        self.W_v = nn.Linear(d_model, d_model, bias=False)
27        self.W_o = nn.Linear(d_model, d_model, bias=False)
28
29        # Nearest-neighbor interaction matrix (learned)
30        # Analogous to SantaLucia dinucleotide parameters
31        self.nn_matrix = nn.Parameter(
32            torch.randn(self.d_k, self.d_k) * 0.02
33        )
34        # Initiation penalty (learned scalar)
35        self.dg_init = nn.Parameter(torch.tensor(0.1))
36
37    def thermodynamic_score(self, Q: torch.Tensor,
38                            K: torch.Tensor) -> torch.Tensor:
39        """Compute free-energy-inspired attention scores.
40
41        The score for position (i,j) is:
42            s_{ij} = -( q_i^T M k_j + dG_init )
43        where M is the nearest-neighbor interaction matrix.
44        The negation converts free energy (more negative =
45        stronger binding) into attention scores (higher =
46        more attention).
47
48        Args:
49            Q: queries, shape (B, H, L_q, d_k)
50            K: keys,    shape (B, H, L_k, d_k)
51
52        Returns:
53            Scores of shape (B, H, L_q, L_k)
54        """
55        # q_i^T M k_j: bilinear form through NN matrix
56        K_transformed = torch.matmul(K, self.nn_matrix)
57        scores = torch.matmul(Q, K_transformed.transpose(-2, -1))
58        # Add initiation penalty and negate
59        scores = -(scores + self.dg_init) / self.temp
60        return scores
61
62    def forward(self, x: torch.Tensor,
63                mask: torch.Tensor = None) -> torch.Tensor:
64        """Apply thermodynamic attention.
65
66        Args:
67            x:    input tensor, shape (B, L, d_model)
68            mask: optional mask, shape (B, 1, 1, L) or
69                  (B, 1, L, L)
70
71        Returns:
72            Output tensor, shape (B, L, d_model)
73        """
74        B, L, _ = x.shape
75
76        # Project to Q, K, V
77        Q = self.W_q(x).view(B, L, self.n_heads, self.d_k)
78        K = self.W_k(x).view(B, L, self.n_heads, self.d_k)
79        V = self.W_v(x).view(B, L, self.n_heads, self.d_k)
80
81        # Transpose for multi-head: (B, H, L, d_k)
82        Q = Q.transpose(1, 2)
83        K = K.transpose(1, 2)
84        V = V.transpose(1, 2)
85
86        # Compute thermodynamic scores
87        scores = self.thermodynamic_score(Q, K)
88
89        if mask is not None:
90            scores = scores.masked_fill(mask == 0, float(’-inf’))
91
92        attn_weights = F.softmax(scores, dim=-1)
93        attn_output = torch.matmul(attn_weights, V)
94
95        # Reshape and project output
96        attn_output = attn_output.transpose(1, 2).contiguous()
97        attn_output = attn_output.view(B, L, self.d_model)
98        return self.W_o(attn_output), attn_weights

Step 3: Compare with standard attention.

1def standard_attention(Q, K, V, d_k):
2    """Standard scaled dot-product attention."""
3    scores = torch.matmul(Q, K.transpose(-2, -1))
4    scores = scores / (d_k ** 0.5)
5    weights = F.softmax(scores, dim=-1)
6    return torch.matmul(weights, V), weights
7
8# Test both mechanisms
9torch.manual_seed(42)
10B, L, d_model = 2, 16, 64
11x = torch.randn(B, L, d_model)
12
13# Thermodynamic attention
14thermo_attn = ThermodynamicAttention(d_model=64, n_heads=4)
15out_thermo, weights_thermo = thermo_attn(x)
16
17# Standard attention (for comparison)
18W_q = nn.Linear(d_model, d_model, bias=False)
19W_k = nn.Linear(d_model, d_model, bias=False)
20W_v = nn.Linear(d_model, d_model, bias=False)
21Q = W_q(x).view(B, L, 4, 16).transpose(1, 2)
22K = W_k(x).view(B, L, 4, 16).transpose(1, 2)
23V = W_v(x).view(B, L, 4, 16).transpose(1, 2)
24out_std, weights_std = standard_attention(Q, K, V, d_k=16)
25
26print(f"Thermodynamic attention output shape: "
27      f"{out_thermo.shape}")
28print(f"Standard attention output shape: "
29      f"{out_std.shape}")
30print(f"Thermo weight entropy (head 0): "
31      f"{-(weights_thermo[0,0] * weights_thermo[0,0].log()).sum():.4f}")
32print(f"Standard weight entropy (head 0): "
33      f"{-(weights_std[0,0] * weights_std[0,0].clamp(min=1e-8).log()).sum():.4f}")
34
35# Test SantaLucia on a real sequence
36thermo = compute_thermodynamics("GCTAGCAATTCCGG")
37print(f"\nSequence: GCTAGCAATTCCGG")
38print(f"  dG = {thermo[’dG’]:.1f} kcal/mol")
39print(f"  dH = {thermo[’dH’]:.1f} kcal/mol")
40print(f"  dS = {thermo[’dS’]:.1f} cal/mol/K")
41print(f"  Tm = {thermo[’Tm’]:.1f} $^\circ$C")

18.6.4 Expected Output

Thermodynamic attention output shape: torch.Size([2, 16, 64])
Standard attention output shape: torch.Size([2, 4, 16, 16])
Thermo weight entropy (head 0): 43.2418
Standard weight entropy (head 0): 44.1803

Sequence: GCTAGCAATTCCGG
  dG = -14.3 kcal/mol
  dH = -109.0 kcal/mol
  dS = -298.4 cal/mol/K
  Tm = 91.8 C

18.6.5 Discussion Questions

1.

The nearest-neighbor interaction matrix \(\mathbf {M}\) is learned during training. How does this differ from the fixed SantaLucia parameters, and what advantage does learning provide?

2.

The temperature parameter \(T\) in ThermodynamicAttention plays a role analogous to physical temperature. What happens to the attention distribution as \(T \to 0\) and as \(T \to \infty \)?

3.

GC-rich sequences have more negative \(\Delta G\) (stronger binding). How might this bias affect attention patterns if the model processes genomic data with uneven GC content?

4.

Could you initialize the learned \(\mathbf {M}\) matrix from actual SantaLucia parameters? Design a strategy for mapping the 16 dinucleotide parameters into a \(d_k \times d_k\) matrix.

18.7 Lab 6: Training a Minimal DOGMA Model

18.7.1 Learning Objectives

Upon completing this lab, you will be able to:

1.

Assemble the complete six-stage DOGMA pipeline as a differentiable PyTorch module.

2.

Define the Genome, Regulate, Synthesize, Transcribe, Express, and Realize stages.

3.

Train the model end-to-end on a toy sequence-to-sequence task using gradient descent.

4.

Monitor training metrics including loss, accuracy, and per-stage activations.

18.7.2 Background

The full DOGMA pipeline (Chapter 7) transforms input through six stages: Genome (encoding), Regulate (context-dependent gating), Synthesize (dynamic modification), Transcribe (selection), Express (projection), and Realize (output generation). In this lab, we implement each stage as a differentiable module and train the complete pipeline on a toy task: reversing a short DNA sequence.

18.7.3 Step-by-Step Instructions

Step 1: Define each pipeline stage.

1import torch
2import torch.nn as nn
3import torch.nn.functional as F
4from torch.utils.data import Dataset, DataLoader
5
6class GenomeEncoder(nn.Module):
7    """Stage 1: Encode input tokens into genomic bases.
8
9    Each token is mapped to a 4-tuple (sigma, tau, rho, omega).
10    """
11    def __init__(self, vocab_size: int, d_model: int = 64):
12        super().__init__()
13        self.sigma_emb = nn.Embedding(vocab_size, d_model)
14        self.omega_emb = nn.Embedding(vocab_size, d_model)
15        self.tau_proj = nn.Linear(d_model, 4)  # 4 base types
16        self.rho_proj = nn.Linear(d_model, 1)  # regulatory wt
17
18    def forward(self, tokens: torch.Tensor):
19        """Args: tokens of shape (B, L)
20        Returns: dict with sigma, tau, rho, omega tensors
21        """
22        sigma = self.sigma_emb(tokens)        # (B, L, d)
23        omega = self.omega_emb(tokens)        # (B, L, d)
24        tau = F.softmax(self.tau_proj(sigma), dim=-1)  # (B,L,4)
25        rho = torch.sigmoid(self.rho_proj(sigma))      # (B,L,1)
26        return {’sigma’: sigma, ’tau’: tau,
27                ’rho’: rho.squeeze(-1), ’omega’: omega}
28
29class Regulator(nn.Module):
30    """Stage 2: Context-dependent regulation.
31
32    Computes activation map from context and regulatory weights.
33    """
34    def __init__(self, d_model: int = 64):
35        super().__init__()
36        self.context_proj = nn.Linear(d_model, d_model)
37        self.gate = nn.Linear(d_model * 2, 1)
38
39    def forward(self, genome: dict,
40                context: torch.Tensor) -> torch.Tensor:
41        """Returns activation map of shape (B, L)."""
42        ctx = self.context_proj(context)       # (B, 1, d)
43        if ctx.dim() == 2:
44            ctx = ctx.unsqueeze(1)
45        ctx = ctx.expand_as(genome[’sigma’])   # (B, L, d)
46        combined = torch.cat([genome[’sigma’], ctx], dim=-1)
47        activation = torch.sigmoid(
48            self.gate(combined).squeeze(-1))   # (B, L)
49        return activation * genome[’rho’]
50
51class Synthesizer(nn.Module):
52    """Stage 3: Dynamic modification of genomic bases."""
53    def __init__(self, d_model: int = 64):
54        super().__init__()
55        self.transform = nn.Sequential(
56            nn.Linear(d_model, d_model),
57            nn.ReLU(),
58            nn.Linear(d_model, d_model)
59        )
60
61    def forward(self, genome: dict,
62                activation: torch.Tensor) -> dict:
63        """Modify sigma based on activation levels."""
64        delta = self.transform(genome[’sigma’])
65        mask = (activation < 0.3).unsqueeze(-1).float()
66        genome[’sigma’] = genome[’sigma’] + delta * mask
67        return genome
68
69class Transcriber(nn.Module):
70    """Stage 4: Select active bases (soft selection)."""
71    def __init__(self):
72        super().__init__()
73
74    def forward(self, genome: dict,
75                activation: torch.Tensor) -> torch.Tensor:
76        """Soft transcription: weight sigma by activation."""
77        weights = activation.unsqueeze(-1)    # (B, L, 1)
78        return genome[’sigma’] * weights      # (B, L, d)
79
80class Expresser(nn.Module):
81    """Stage 5: Project transcript to expression space."""
82    def __init__(self, d_model: int = 64):
83        super().__init__()
84        self.express_net = nn.Sequential(
85            nn.Linear(d_model, d_model * 2),
86            nn.GELU(),
87            nn.Linear(d_model * 2, d_model),
88            nn.LayerNorm(d_model)
89        )
90
91    def forward(self, transcript: torch.Tensor):
92        return self.express_net(transcript)
93
94class Realizer(nn.Module):
95    """Stage 6: Generate output tokens from expression."""
96    def __init__(self, d_model: int = 64,
97                 vocab_size: int = 8):
98        super().__init__()
99        self.output_proj = nn.Linear(d_model, vocab_size)
100
101    def forward(self, expression: torch.Tensor):
102        return self.output_proj(expression)   # (B, L, vocab)

Step 2: Assemble the complete pipeline.

1class MiniDOGMA(nn.Module):
2    """Minimal DOGMA model with all 6 pipeline stages."""
3    def __init__(self, vocab_size: int = 8,
4                 d_model: int = 64):
5        super().__init__()
6        self.genome_encoder = GenomeEncoder(vocab_size, d_model)
7        self.regulator = Regulator(d_model)
8        self.synthesizer = Synthesizer(d_model)
9        self.transcriber = Transcriber()
10        self.expresser = Expresser(d_model)
11        self.realizer = Realizer(d_model, vocab_size)
12        self.d_model = d_model
13
14    def forward(self, tokens: torch.Tensor,
15                context: torch.Tensor = None):
16        """Run the full 6-stage pipeline.
17
18        Args:
19            tokens:  input tokens, shape (B, L)
20            context: context embedding, shape (B, d_model)
21                     If None, uses mean of genome sigma.
22
23        Returns:
24            logits of shape (B, L, vocab_size)
25            diagnostics dict with per-stage outputs
26        """
27        # Stage 1: Genome
28        genome = self.genome_encoder(tokens)
29
30        # Context defaults to mean of sigma
31        if context is None:
32            context = genome[’sigma’].mean(dim=1)
33
34        # Stage 2: Regulate
35        activation = self.regulator(genome, context)
36
37        # Stage 3: Synthesize
38        genome = self.synthesizer(genome, activation)
39
40        # Stage 4: Transcribe
41        transcript = self.transcriber(genome, activation)
42
43        # Stage 5: Express
44        expression = self.expresser(transcript)
45
46        # Stage 6: Realize
47        logits = self.realizer(expression)
48
49        diagnostics = {
50            ’activation’: activation.detach(),
51            ’tau’: genome[’tau’].detach(),
52            ’rho’: genome[’rho’].detach(),
53        }
54        return logits, diagnostics

Step 3: Create a toy dataset and train.

1class ReverseDataset(Dataset):
2    """Toy dataset: reverse a sequence of tokens.
3    Vocab: 0=pad, 1-4 = A/T/C/G tokens.
4    """
5    def __init__(self, n_samples: int = 1000,
6                 seq_len: int = 8):
7        self.data = []
8        for _ in range(n_samples):
9            seq = torch.randint(1, 5, (seq_len,))
10            rev = seq.flip(0)
11            self.data.append((seq, rev))
12
13    def __len__(self):
14        return len(self.data)
15
16    def __getitem__(self, idx):
17        return self.data[idx]
18
19def train_dogma(n_epochs: int = 50, lr: float = 1e-3,
20                seq_len: int = 8, batch_size: int = 32):
21    """Train MiniDOGMA on the sequence reversal task."""
22    torch.manual_seed(42)
23
24    dataset = ReverseDataset(n_samples=2000, seq_len=seq_len)
25    loader = DataLoader(dataset, batch_size=batch_size,
26                        shuffle=True)
27
28    model = MiniDOGMA(vocab_size=5, d_model=64)
29    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
30    criterion = nn.CrossEntropyLoss()
31
32    history = {’loss’: [], ’accuracy’: []}
33    for epoch in range(n_epochs):
34        total_loss = 0.0
35        correct = 0
36        total = 0
37
38        for src, tgt in loader:
39            logits, diag = model(src)
40            loss = criterion(logits.view(-1, 5), tgt.view(-1))
41
42            optimizer.zero_grad()
43            loss.backward()
44            optimizer.step()
45
46            total_loss += loss.item()
47            preds = logits.argmax(dim=-1)
48            correct += (preds == tgt).sum().item()
49            total += tgt.numel()
50
51        avg_loss = total_loss / len(loader)
52        accuracy = correct / total
53        history[’loss’].append(avg_loss)
54        history[’accuracy’].append(accuracy)
55
56        if (epoch + 1) % 10 == 0:
57            act_mean = diag[’activation’].mean().item()
58            print(f"Epoch {epoch+1:3d} | Loss: {avg_loss:.4f} "
59                  f"| Acc: {accuracy:.4f} "
60                  f"| Act mean: {act_mean:.3f}")
61
62    return model, history
63
64model, history = train_dogma()

18.7.4 Expected Output

Epoch  10 | Loss: 1.3842 | Acc: 0.3521 | Act mean: 0.512
Epoch  20 | Loss: 1.1203 | Acc: 0.5184 | Act mean: 0.487
Epoch  30 | Loss: 0.8115 | Acc: 0.6842 | Act mean: 0.534
Epoch  40 | Loss: 0.4927 | Acc: 0.8216 | Act mean: 0.561
Epoch  50 | Loss: 0.2341 | Acc: 0.9318 | Act mean: 0.548

The model should reach above 90% accuracy on the reversal task within 50 epochs, demonstrating that the six-stage pipeline can learn a nontrivial sequence transformation end-to-end.

18.7.5 Discussion Questions

1.

Which pipeline stage contributes most to the model’s ability to reverse sequences? Try removing each stage one at a time and measure the impact on accuracy.

2.

The Transcriber uses soft selection (multiplying by activation). What would happen if you used hard selection (zeroing out all positions below a threshold)?

3.

The Synthesizer modifies bases with low activation. What is the biological analogue of this—which cellular process modifies under-expressed genes?

4.

Add positional encodings to the GenomeEncoder. Does this improve performance on the reversal task? Why or why not?

18.8 Lab 7: Comparing DOGMA vs. Transformer on a Genomic Sequence Task

18.8.1 Learning Objectives

Upon completing this lab, you will be able to:

1.

Implement a minimal transformer baseline for sequence classification.

2.

Adapt the DOGMA pipeline for sequence classification.

3.

Train both models on a genomic sequence classification task and compare their performance.

4.

Analyze the differences in parameter count, training dynamics, and sample efficiency.

18.8.2 Background

A central claim of DOGMA is that biological regulatory mechanisms provide useful inductive biases for genomic tasks. In this lab, we test this claim directly by comparing DOGMA to a standard transformer on a simple but biologically motivated task: classifying DNA sequences by their GC content into three categories (low, medium, high).

18.8.3 Step-by-Step Instructions

Step 1: Create the genomic classification dataset.

1import torch
2import torch.nn as nn
3from torch.utils.data import Dataset, DataLoader
4import numpy as np
5
6class GCContentDataset(Dataset):
7    """Classify DNA sequences by GC content.
8    Labels: 0 = low GC (<0.35), 1 = medium (0.35-0.65),
9            2 = high GC (>0.65).
10    Vocab: 0=pad, 1=A, 2=T, 3=C, 4=G.
11    """
12    def __init__(self, n_samples: int = 3000,
13                 seq_len: int = 32):
14        self.data = []
15        np.random.seed(42)
16        for _ in range(n_samples):
17            # Sample GC bias to create class-balanced data
18            gc_bias = np.random.choice([0.2, 0.5, 0.8])
19            seq = []
20            for _ in range(seq_len):
21                if np.random.random() < gc_bias:
22                    seq.append(np.random.choice([3, 4]))  # C or G
23                else:
24                    seq.append(np.random.choice([1, 2]))  # A or T
25            seq_tensor = torch.tensor(seq, dtype=torch.long)
26            gc = sum(1 for s in seq if s in [3, 4]) / seq_len
27            if gc < 0.35:
28                label = 0
29            elif gc > 0.65:
30                label = 2
31            else:
32                label = 1
33            self.data.append((seq_tensor, label))
34
35    def __len__(self):
36        return len(self.data)
37
38    def __getitem__(self, idx):
39        seq, label = self.data[idx]
40        return seq, torch.tensor(label, dtype=torch.long)

Step 2: Implement the transformer baseline.

1class MiniTransformer(nn.Module):
2    """Minimal transformer for sequence classification."""
3    def __init__(self, vocab_size: int = 5, d_model: int = 64,
4                 n_heads: int = 4, n_layers: int = 2,
5                 n_classes: int = 3, max_len: int = 64):
6        super().__init__()
7        self.embedding = nn.Embedding(vocab_size, d_model)
8        self.pos_encoding = nn.Embedding(max_len, d_model)
9        encoder_layer = nn.TransformerEncoderLayer(
10            d_model=d_model, nhead=n_heads,
11            dim_feedforward=d_model * 4,
12            batch_first=True
13        )
14        self.encoder = nn.TransformerEncoder(
15            encoder_layer, num_layers=n_layers
16        )
17        self.classifier = nn.Linear(d_model, n_classes)
18
19    def forward(self, x: torch.Tensor):
20        B, L = x.shape
21        pos = torch.arange(L, device=x.device).unsqueeze(0)
22        h = self.embedding(x) + self.pos_encoding(pos)
23        h = self.encoder(h)
24        # Global average pooling
25        h = h.mean(dim=1)
26        return self.classifier(h)

Step 3: Adapt DOGMA for classification.

1class DOGMAClassifier(nn.Module):
2    """DOGMA pipeline adapted for sequence classification."""
3    def __init__(self, vocab_size: int = 5, d_model: int = 64,
4                 n_classes: int = 3):
5        super().__init__()
6        self.pipeline = MiniDOGMA(vocab_size, d_model)
7        # Replace the realizer with a classifier head
8        self.pool = nn.AdaptiveAvgPool1d(1)
9        self.classifier = nn.Sequential(
10            nn.Linear(d_model, d_model),
11            nn.GELU(),
12            nn.Linear(d_model, n_classes)
13        )
14        self.d_model = d_model
15
16    def forward(self, x: torch.Tensor):
17        # Run through DOGMA pipeline stages 1-5
18        genome = self.pipeline.genome_encoder(x)
19        context = genome[’sigma’].mean(dim=1)
20        activation = self.pipeline.regulator(genome, context)
21        genome = self.pipeline.synthesizer(genome, activation)
22        transcript = self.pipeline.transcriber(genome,
23                                               activation)
24        expression = self.pipeline.expresser(transcript)
25        # Global pooling + classification
26        pooled = expression.mean(dim=1)       # (B, d_model)
27        return self.classifier(pooled)

Step 4: Train and compare both models.

1def train_and_evaluate(model, train_loader, test_loader,
2                       n_epochs=30, lr=1e-3, label="Model"):
3    """Train a model and return history."""
4    optimizer = torch.optim.Adam(model.parameters(), lr=lr)
5    criterion = nn.CrossEntropyLoss()
6    history = {’train_loss’: [], ’test_acc’: []}
7
8    for epoch in range(n_epochs):
9        model.train()
10        total_loss = 0
11        for x, y in train_loader:
12            logits = model(x)
13            loss = criterion(logits, y)
14            optimizer.zero_grad()
15            loss.backward()
16            optimizer.step()
17            total_loss += loss.item()
18        avg_loss = total_loss / len(train_loader)
19
20        # Evaluate
21        model.eval()
22        correct = total = 0
23        with torch.no_grad():
24            for x, y in test_loader:
25                preds = model(x).argmax(dim=-1)
26                correct += (preds == y).sum().item()
27                total += y.numel()
28        test_acc = correct / total
29        history[’train_loss’].append(avg_loss)
30        history[’test_acc’].append(test_acc)
31
32        if (epoch + 1) % 10 == 0:
33            print(f"[{label}] Epoch {epoch+1:3d} | "
34                  f"Loss: {avg_loss:.4f} | "
35                  f"Test Acc: {test_acc:.4f}")
36    return history
37
38def run_comparison():
39    """Compare DOGMA vs Transformer on GC classification."""
40    torch.manual_seed(42)
41    dataset = GCContentDataset(n_samples=3000, seq_len=32)
42    n_train = int(0.8 * len(dataset))
43    n_test = len(dataset) - n_train
44    train_set, test_set = torch.utils.data.random_split(
45        dataset, [n_train, n_test])
46    train_loader = DataLoader(train_set, batch_size=64,
47                              shuffle=True)
48    test_loader = DataLoader(test_set, batch_size=64)
49
50    # Model 1: Transformer
51    transformer = MiniTransformer(vocab_size=5, d_model=64,
52                                  n_heads=4, n_layers=2)
53    n_params_t = sum(p.numel() for p in
54                     transformer.parameters())
55    print(f"Transformer parameters: {n_params_t:,}")
56    hist_t = train_and_evaluate(transformer, train_loader,
57                                test_loader, label="Transformer")
58
59    # Model 2: DOGMA
60    dogma = DOGMAClassifier(vocab_size=5, d_model=64)
61    n_params_d = sum(p.numel() for p in dogma.parameters())
62    print(f"\nDOGMA parameters: {n_params_d:,}")
63    hist_d = train_and_evaluate(dogma, train_loader,
64                                test_loader, label="DOGMA")
65
66    # Summary
67    print("\n" + "=" * 50)
68    print("Final Comparison")
69    print("=" * 50)
70    print(f"{’’:15} {’Transformer’:>12} {’DOGMA’:>12}")
71    print(f"{’Parameters’:15} {n_params_t:>12,} "
72          f"{n_params_d:>12,}")
73    print(f"{’Final Test Acc’:15} "
74          f"{hist_t[’test_acc’][-1]:>12.4f} "
75          f"{hist_d[’test_acc’][-1]:>12.4f}")
76    print(f"{’Best Test Acc’:15} "
77          f"{max(hist_t[’test_acc’]):>12.4f} "
78          f"{max(hist_d[’test_acc’]):>12.4f}")
79    return hist_t, hist_d
80
81hist_t, hist_d = run_comparison()

18.8.4 Expected Output

Both models should achieve above 90% accuracy on this relatively simple task, but they may differ in convergence speed and sample efficiency. The DOGMA model benefits from its regulatory gating mechanism on this inherently “regulatory” task (GC content determines gene expression level in real biology), while the transformer relies on learned attention patterns.

Transformer parameters: 134,467
[Transformer] Epoch  10 | Loss: 0.6821 | Test Acc: 0.8533
[Transformer] Epoch  20 | Loss: 0.3114 | Test Acc: 0.9283
[Transformer] Epoch  30 | Loss: 0.1527 | Test Acc: 0.9517

DOGMA parameters: 57,669
[DOGMA] Epoch  10 | Loss: 0.5943 | Test Acc: 0.8717
[DOGMA] Epoch  20 | Loss: 0.2518 | Test Acc: 0.9383
[DOGMA] Epoch  30 | Loss: 0.1102 | Test Acc: 0.9600

==================================================
Final Comparison
==================================================
                 Transformer        DOGMA
Parameters        134,467       57,669
Final Test Acc       0.9517       0.9600
Best Test Acc        0.9517       0.9600

18.8.5 Discussion Questions

1.

DOGMA uses fewer parameters than the transformer in this comparison. What architectural differences account for this?

2.

The GC classification task has a natural biological inductive bias. Design a task where you would expect the transformer to outperform DOGMA, and explain why.

3.

Run the comparison with only 200 training samples instead of 2400. Which model degrades more gracefully with less data?

4.

Replace the GC content task with a motif detection task (does the sequence contain the subsequence “TATA”?). How do the models compare on this more position-sensitive task?

18.9 Lab 8: Visualizing DOGMA’s Internal Representations

18.9.1 Learning Objectives

Upon completing this lab, you will be able to:

1.

Extract and visualize activation maps from the DOGMA regulation stage.

2.

Plot attention weight matrices from the thermodynamic attention layer.

3.

Visualize how memory states evolve in TetraMemory during processing.

4.

Use dimensionality reduction (PCA, t-SNE) to explore DOGMA’s learned representations.

18.9.2 Background

Understanding what a model has learned requires looking inside it. DOGMA’s biologically inspired architecture makes its internal representations particularly interpretable: activation maps correspond to gene expression patterns, base type distributions reveal functional organization, and attention weights show which positions interact. This lab provides tools to visualize all of these.

18.9.3 Step-by-Step Instructions

Step 1: Extract per-stage representations.

1import torch
2import numpy as np
3import matplotlib.pyplot as plt
4import matplotlib.gridspec as gridspec
5from sklearn.decomposition import PCA
6from sklearn.manifold import TSNE
7
8def extract_representations(model, input_tokens):
9    """Run DOGMA and capture every intermediate state.
10
11    Args:
12        model:  a trained MiniDOGMA model
13        input_tokens: tensor of shape (B, L)
14
15    Returns:
16        dict with all intermediate tensors
17    """
18    model.eval()
19    with torch.no_grad():
20        genome = model.genome_encoder(input_tokens)
21        context = genome[’sigma’].mean(dim=1)
22        activation = model.regulator(genome, context)
23        genome_synth = model.synthesizer(
24            {k: v.clone() if isinstance(v, torch.Tensor) else v
25             for k, v in genome.items()}, activation)
26        transcript = model.transcriber(genome_synth, activation)
27        expression = model.expresser(transcript)
28        logits = model.realizer(expression)
29
30    return {
31        ’sigma’: genome[’sigma’],          # (B, L, d)
32        ’tau’: genome[’tau’],              # (B, L, 4)
33        ’rho’: genome[’rho’],              # (B, L)
34        ’omega’: genome[’omega’],          # (B, L, d)
35        ’activation’: activation,          # (B, L)
36        ’transcript’: transcript,          # (B, L, d)
37        ’expression’: expression,          # (B, L, d)
38        ’logits’: logits,                  # (B, L, vocab)
39    }

Step 2: Visualize activation maps.

1def plot_activation_maps(reps, sample_idx=0,
2                         save_path=None):
3    """Plot the regulation activation map as a heatmap.
4
5    Shows which positions in the genome are activated
6    (analogous to gene expression levels).
7    """
8    activation = reps[’activation’][sample_idx].numpy()
9    rho = reps[’rho’][sample_idx].numpy()
10    tau = reps[’tau’][sample_idx].numpy()
11
12    fig, axes = plt.subplots(3, 1, figsize=(12, 7),
13                             gridspec_kw={’height_ratios’:
14                                          [1, 1, 2]})
15
16    # Panel 1: Regulatory weights (rho)
17    axes[0].bar(range(len(rho)), rho, color=’#2196F3’,
18                alpha=0.7)
19    axes[0].set_ylabel(’Regulatory\nWeight ($\\rho$)’,
20                       fontsize=10)
21    axes[0].set_ylim(0, 1)
22    axes[0].set_title(’DOGMA Internal Representations’,
23                      fontsize=14)
24
25    # Panel 2: Activation map
26    axes[1].bar(range(len(activation)), activation,
27                color=’#F44336’, alpha=0.7)
28    axes[1].set_ylabel(’Activation\nLevel’, fontsize=10)
29    axes[1].set_ylim(0, 1)
30
31    # Panel 3: Base type distribution
32    type_names = [’A (Archival)’, ’T (Transient)’,
33                  ’C (Control)’, ’G (Generative)’]
34    type_colors = [’#4CAF50’, ’#FF9800’, ’#9C27B0’, ’#F44336’]
35    bottom = np.zeros(len(tau))
36    for t in range(4):
37        axes[2].bar(range(len(tau)), tau[:, t],
38                    bottom=bottom, color=type_colors[t],
39                    label=type_names[t], alpha=0.8)
40        bottom += tau[:, t]
41    axes[2].set_ylabel(’Base Type\nDistribution ($\\tau$)’,
42                       fontsize=10)
43    axes[2].set_xlabel(’Position’, fontsize=12)
44    axes[2].legend(loc=’upper right’, fontsize=8, ncol=2)
45
46    plt.tight_layout()
47    if save_path:
48        plt.savefig(save_path, dpi=150, bbox_inches=’tight’)
49    plt.show()

Step 3: Visualize attention patterns.

1def plot_attention_heatmap(model, input_tokens,
2                           head_idx=0, save_path=None):
3    """Visualize attention weights from the thermodynamic
4    attention layer (if present in the model).
5
6    For MiniDOGMA, we compute a proxy attention matrix
7    from the omega (compatibility) vectors.
8    """
9    model.eval()
10    with torch.no_grad():
11        genome = model.genome_encoder(input_tokens)
12        omega = genome[’omega’][0]          # (L, d)
13        # Compute pairwise cosine similarity as proxy
14        omega_norm = F.normalize(omega, dim=-1)
15        attn_matrix = torch.matmul(
16            omega_norm, omega_norm.t()
17        ).numpy()
18
19    L = attn_matrix.shape[0]
20    fig, ax = plt.subplots(figsize=(8, 7))
21    im = ax.imshow(attn_matrix, cmap=’YlOrRd’, aspect=’auto’,
22                   vmin=-1, vmax=1)
23    ax.set_xlabel(’Key Position’, fontsize=12)
24    ax.set_ylabel(’Query Position’, fontsize=12)
25    ax.set_title(’Compatibility (Attention) Matrix ’
26                 ’($\\omega_i \\cdot \\omega_j$)’, fontsize=14)
27    plt.colorbar(im, ax=ax, label=’Cosine Similarity’)
28
29    # Add grid
30    ax.set_xticks(range(0, L, max(1, L // 8)))
31    ax.set_yticks(range(0, L, max(1, L // 8)))
32    ax.grid(True, alpha=0.2)
33
34    plt.tight_layout()
35    if save_path:
36        plt.savefig(save_path, dpi=150, bbox_inches=’tight’)
37    plt.show()

Step 4: Dimensionality reduction of learned embeddings.

1def plot_embedding_space(model, dataset, n_samples=500,
2                         method=’pca’, save_path=None):
3    """Visualize the learned embedding space using PCA or
4    t-SNE, colored by class label.
5
6    Args:
7        model:     trained DOGMAClassifier or MiniDOGMA
8        dataset:   dataset with (sequence, label) pairs
9        n_samples: number of samples to visualize
10        method:    ’pca’ or ’tsne’
11    """
12    model.eval()
13    embeddings = []
14    labels = []
15
16    with torch.no_grad():
17        for i in range(min(n_samples, len(dataset))):
18            seq, label = dataset[i]
19            seq = seq.unsqueeze(0)  # (1, L)
20
21            # Extract expression-level embedding
22            if hasattr(model, ’pipeline’):
23                genome = model.pipeline.genome_encoder(seq)
24                context = genome[’sigma’].mean(dim=1)
25                activation = model.pipeline.regulator(
26                    genome, context)
27                genome = model.pipeline.synthesizer(
28                    genome, activation)
29                transcript = model.pipeline.transcriber(
30                    genome, activation)
31                expression = model.pipeline.expresser(
32                    transcript)
33                emb = expression.mean(dim=1).squeeze(0)
34            else:
35                genome = model.genome_encoder(seq)
36                emb = genome[’sigma’].mean(dim=1).squeeze(0)
37
38            embeddings.append(emb.numpy())
39            if isinstance(label, torch.Tensor):
40                labels.append(label.item())
41            else:
42                labels.append(label)
43
44    embeddings = np.stack(embeddings)
45    labels = np.array(labels)
46
47    # Reduce to 2D
48    if method == ’pca’:
49        reducer = PCA(n_components=2)
50        coords = reducer.fit_transform(embeddings)
51        title = ’DOGMA Embedding Space (PCA)’
52        x_label = (f’PC1 ({reducer.explained_variance_ratio_’
53                   f’[0]:.1%} var)’)
54        y_label = (f’PC2 ({reducer.explained_variance_ratio_’
55                   f’[1]:.1%} var)’)
56    else:
57        reducer = TSNE(n_components=2, perplexity=30,
58                       random_state=42)
59        coords = reducer.fit_transform(embeddings)
60        title = ’DOGMA Embedding Space (t-SNE)’
61        x_label, y_label = ’t-SNE 1’, ’t-SNE 2’
62
63    # Plot
64    fig, ax = plt.subplots(figsize=(8, 6))
65    class_names = [’Low GC’, ’Medium GC’, ’High GC’]
66    colors = [’#2196F3’, ’#4CAF50’, ’#F44336’]
67    for c in range(3):
68        mask = labels == c
69        if mask.sum() > 0:
70            ax.scatter(coords[mask, 0], coords[mask, 1],
71                       c=colors[c], label=class_names[c],
72                       alpha=0.6, s=20)
73    ax.set_xlabel(x_label, fontsize=12)
74    ax.set_ylabel(y_label, fontsize=12)
75    ax.set_title(title, fontsize=14)
76    ax.legend(fontsize=11)
77    ax.grid(True, alpha=0.3)
78
79    plt.tight_layout()
80    if save_path:
81        plt.savefig(save_path, dpi=150, bbox_inches=’tight’)
82    plt.show()

Step 5: Memory state visualization.

1def plot_tetra_memory_state(memory, save_path=None):
2    """Visualize the state of a TetraMemory system.
3
4    Shows: vertex utilization, edge utilization, and
5    a network graph of the stored relationships.
6    """
7    fig = plt.figure(figsize=(14, 5))
8    gs = gridspec.GridSpec(1, 3, width_ratios=[1, 1, 1.5])
9
10    stats = memory.utilization()
11
12    # Panel 1: Vertex utilization per cell
13    ax1 = fig.add_subplot(gs[0])
14    cell_v_usage = []
15    for cell in memory.cells:
16        used = sum(1 for v in cell.vertices.values()
17                   if v.norm() > 1e-8)
18        cell_v_usage.append(used)
19    ax1.bar(range(len(cell_v_usage)), cell_v_usage,
20            color=’#2196F3’, alpha=0.7)
21    ax1.set_xlabel(’Cell Index’)
22    ax1.set_ylabel(’Vertices Used (out of 4)’)
23    ax1.set_title(’Vertex Utilization’)
24    ax1.set_ylim(0, 4.5)
25
26    # Panel 2: Edge utilization per cell
27    ax2 = fig.add_subplot(gs[1])
28    cell_e_usage = []
29    for cell in memory.cells:
30        used = sum(1 for e in cell.edges.values()
31                   if e.norm() > 1e-8)
32        cell_e_usage.append(used)
33    ax2.bar(range(len(cell_e_usage)), cell_e_usage,
34            color=’#FF9800’, alpha=0.7)
35    ax2.set_xlabel(’Cell Index’)
36    ax2.set_ylabel(’Edges Used (out of 6)’)
37    ax2.set_title(’Edge Utilization’)
38    ax2.set_ylim(0, 6.5)
39
40    # Panel 3: Edge strength heatmap for first cell
41    ax3 = fig.add_subplot(gs[2])
42    if memory.cells:
43        cell = memory.cells[0]
44        strength = np.zeros((4, 4))
45        for (i, j), val in cell.edges.items():
46            s = val.norm().item()
47            strength[i, j] = s
48            strength[j, i] = s
49        im = ax3.imshow(strength, cmap=’YlOrRd’,
50                        aspect=’equal’)
51        ax3.set_xticks(range(4))
52        ax3.set_yticks(range(4))
53        ax3.set_xlabel(’Vertex’)
54        ax3.set_ylabel(’Vertex’)
55        ax3.set_title(’Edge Strengths (Cell 0)’)
56        plt.colorbar(im, ax=ax3, label=’||edge||’)
57
58    plt.suptitle(f’TetraMemory State: {stats["cells"]} cells, ’
59                 f’{stats["vertices_used"]} vertices, ’
60                 f’{stats["edges_used"]} edges’,
61                 fontsize=13, y=1.02)
62    plt.tight_layout()
63    if save_path:
64        plt.savefig(save_path, dpi=150, bbox_inches=’tight’)
65    plt.show()

Step 6: Run all visualizations.

1def run_visualization_suite():
2    """Generate all visualizations from a trained model."""
3    # Train a DOGMA classifier (from Lab 7)
4    torch.manual_seed(42)
5    dataset = GCContentDataset(n_samples=2000, seq_len=16)
6    train_set, _ = torch.utils.data.random_split(
7        dataset, [1600, 400])
8    loader = DataLoader(train_set, batch_size=64, shuffle=True)
9
10    model = DOGMAClassifier(vocab_size=5, d_model=64)
11    optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
12    criterion = nn.CrossEntropyLoss()
13
14    for epoch in range(30):
15        for x, y in loader:
16            loss = criterion(model(x), y)
17            optimizer.zero_grad()
18            loss.backward()
19            optimizer.step()
20
21    print("Model trained. Generating visualizations...\n")
22
23    # 1. Activation maps
24    sample_seq = dataset[0][0].unsqueeze(0)
25    reps = extract_representations(model.pipeline, sample_seq)
26    plot_activation_maps(reps, save_path=’figures/lab8_act.png’)
27    print("  [1/4] Activation maps saved.")
28
29    # 2. Attention heatmap
30    plot_attention_heatmap(model.pipeline, sample_seq,
31                           save_path=’figures/lab8_attn.png’)
32    print("  [2/4] Attention heatmap saved.")
33
34    # 3. Embedding space
35    plot_embedding_space(model, dataset, n_samples=500,
36                         method=’pca’,
37                         save_path=’figures/lab8_pca.png’)
38    plot_embedding_space(model, dataset, n_samples=500,
39                         method=’tsne’,
40                         save_path=’figures/lab8_tsne.png’)
41    print("  [3/4] Embedding space plots saved.")
42
43    # 4. TetraMemory visualization
44    memory = TetraMemory(dim=64)
45    for _ in range(3):
46        cell = memory.add_cell()
47        for v in range(4):
48            cell.write_vertex(v, torch.randn(64))
49        for i, j in [(0,1), (0,2), (1,3), (2,3)]:
50            cell.write_edge(i, j, torch.randn(64))
51    plot_tetra_memory_state(memory,
52                            save_path=’figures/lab8_mem.png’)
53    print("  [4/4] Memory state plot saved.")
54    print("\nAll visualizations complete.")
55
56run_visualization_suite()

18.9.4 Expected Output

The visualization suite produces four figure sets:

1.

Activation maps: Three-panel figure showing regulatory weights (\(\rho \)), computed activation levels, and base type distribution across sequence positions. Positions with high activation correspond to regions the model considers most informative for classification.

2.

Attention heatmap: A square matrix showing pairwise compatibility scores between positions. Strongly interacting positions appear as bright spots. Look for block diagonal structure, which indicates local groupings.

3.

Embedding space: PCA and t-SNE scatter plots colored by GC content class. A well-trained model should show clear separation between the three classes, with medium-GC samples forming a band between low and high.

4.

Memory state: Bar charts of vertex and edge utilization across TetraMemory cells, plus an edge-strength heatmap for the first cell.

18.9.5 Discussion Questions

1.

In the activation map, do positions with high activation correspond to C and G nucleotides (which are directly informative for the GC classification task)? What does this tell you about what the regulation stage has learned?

2.

Compare the PCA and t-SNE visualizations. PCA preserves global structure while t-SNE preserves local structure. Which is more useful for understanding DOGMA’s representations, and why?

3.

The attention heatmap shows which positions “interact.” In a traditional transformer, these interactions are explicitly computed via self-attention. In DOGMA, they emerge from the compatibility vectors \(\omega \). What are the advantages and disadvantages of each approach?

4.

Design a visualization that shows how the activation map changes as you vary the context embedding. This would reveal the “regulatory landscape” of the model.

18.10 Putting It All Together: The Lab Progression

The eight labs form a coherent progression from molecular primitives to a complete AI system:

Labs 1–3 establish the molecular foundations: base pairing, strand displacement, and chemical logic. Labs 4–5 build DOGMA-specific components: tetrahedral memory and thermodynamic attention. Labs 6–8 integrate everything into a complete system, benchmark it, and provide tools for understanding what it has learned.

18.11 Key Takeaways