Eigenia Digital Twin Complete Mathematical Ontology (40 Governing Equations)

J. McKenney

This is a standalone technical reference for the Cyber Digital Twin engine rather than a numbered entry in a series: it documents the forty mathematical formulas implemented across six computational engines, Monte Carlo walk, edge weight composition, Hawkes cascade, structural causal model, ALE insurance and ATQ actor scoring, that the corpus's Monte Carlo and insurance-quantification papers draw on.

Licence: CC BY 4.0. 17 September 2026.

Executive Abstract#

The Cyber Digital Twin engine runs on forty mathematical formulas spread across six pieces of software: a Monte Carlo random-walk simulator that performs Boltzmann-weighted walks on attack graphs, an edge-weighting calculation that folds a base weight for each relationship type together with fourteen modifier dimensions, a self-exciting event model for correlated losses, an eight-stage cause-and-effect model, an insurance loss-estimation engine, and an actor-scoring procedure. Every formula already lives at a specific file and line of running code.

The eight stages form a consequence cascade, written CS0 through CS7, that traces how a disturbance moves outward from physical process to cyber detection, isolation, organizational and geopolitical effect, economic quantification, psychographic shift, and temporal forecast. This numbering is distinct from the eight-layer digital twin graph and from the completeness levels used elsewhere in the corpus; they share no meaning at any index.

This document exists so the mathematics of a formula and its implementation can be checked against each other. For each one it names the file and line, states the equation in the notation the code uses, and gives the plain meaning of every term.

It makes no argument and defends no result. It is a reference, not a treatise, and its only claim is that the mathematics shown matches the code cited beside it. A reader who wants the reasoning behind a formula rather than its statement is pointed to the paper that uses it.

Abstract#

The Cyber Digital Twin implements forty mathematical formulas across six computational engines. Together they span the full eight-stage consequence cascade, CS0 Physical Process through CS7 Temporal Forecast, connecting physical process dynamics to economic quantification and predictive forecasting. Named constructs include Boltzmann edge selection, Pareto and generalized Pareto tail sampling, the Hill and percentile tail-index estimators, conditional value at risk, the Hawkes intensity function and its reproduction number, Pearl do-calculus stage mechanisms, the Poisson-Pareto annual loss expectancy model with Gordon-Loeb optimal investment, and the eight-component actor threat quotient. Each engine, mc-engine.ts, mc-weights.ts, mc-hawkes.ts, mc-scm.ts, ale-engine.ts and atq-migration.sql, is cited by file and line. All formulas are extracted directly from the production codebase; none derives from external research.

1. Notation and Engine Overview#

The engines and their roles:

All formulas below are extracted directly from the production codebase.

2. Monte Carlo Walk Engine#

F1. Mulberry32 PRNG#

Seedable 32-bit pseudorandom number generator for reproducible simulations.

File: mc-engine.ts:27-36

export function createPRNG(seed?: number | null): () => number {
    if (seed == null) return Math.random;
    let s = seed | 0;
    return () => {
        s |= 0; s = s + 0x6D2B79F5 | 0;
        let t = Math.imul(s ^ s >>> 15, 1 | s);
        t = t + Math.imul(t ^ t >>> 7, 61 | t) ^ t;
        return ((t ^ t >>> 14) >>> 0) / 4294967296;
    };
}

Output range: [0, 1). When seed is null, falls back to Math.random().

F2. Boltzmann Distribution (Edge Selection)#

Selects the next edge during a random walk using a Boltzmann (softmax) distribution with temperature parameter T.

File: mc-engine.ts:320-344

Mathematical notation:

Where w_i is the edge weight and T is the temperature parameter. Higher temperature = more exploration (uniform). Lower temperature = more exploitation (greedy).

// Boltzmann: P(e) = exp(weight / T)
const maxW = Math.max(...candidates.map(e => e.weight));
const energies = candidates.map(e => Math.exp((e.weight - maxW) / temperature)); // log-sum-exp trick
const Z = energies.reduce((a, b) => a + b, 0);
const probs = energies.map(e => e / Z);

The log-sum-exp trick (subtracting maxW) prevents numerical overflow. Black swan walks use T_bs = T * 4.0 for maximum exploration.

F3. Pareto Sampling (Fat-Tail Cost)#

Samples breach costs from a Pareto distribution for fat-tailed loss modeling.

File: mc-engine.ts:348-351

Mathematical notation:

Where U ~ Uniform(0,1), alpha is the tail index, and x_min is the minimum loss threshold.

function paretoSample(alpha: number, xmin: number = 100000, rng: () => number = Math.random): number {
    const u = rng();
    return xmin * Math.pow(1 - u, -1.0 / alpha);
}

F4. Hill Estimator (Tail Index)#

Estimates the Pareto tail index alpha from observed cost data using the top 10 percent order statistics.

File: mc-engine.ts:355-366

Mathematical notation:

Where X_(1) >= X_(2) >= ... >= X_(k) are the top k order statistics and k = floor(0.1 * n).

function hillEstimator(data: number[]): number {
    if (data.length < 10) return 1.5; // default
    const sorted = [...data].sort((a, b) => b - a);
    const k = Math.max(5, Math.floor(data.length * 0.1)); // top 10%
    const xk = sorted[k];
    if (xk <= 0) return 1.5;
    let sum = 0;
    for (let i = 0; i < k; i++) {
        if (sorted[i] > 0 && xk > 0) sum += Math.log(sorted[i] / xk);
    }
    return sum > 0 ? k / sum : 1.5;
}

F5. CVaR (Conditional Value at Risk)#

The expected loss given that the loss exceeds the VaR threshold.

File: mc-engine.ts:705-708

Mathematical notation:

Where T_alpha = {l : l >= VaR_alpha} is the tail set.

const tail95 = allCosts.filter(c => c >= var95);
const cvar95 = tail95.length ? tail95.reduce((a, b) => a + b, 0) / tail95.length : 0;
const tail99 = allCosts.filter(c => c >= var99);
const cvar99 = tail99.length ? tail99.reduce((a, b) => a + b, 0) / tail99.length : 0;

F6. Gaussian vs Pareto Ratio#

Measures how wrong a Gaussian assumption would be. Values > 2.0 indicate the Gaussian model is dangerously inadequate.

File: mc-engine.ts:776-778

Mathematical notation:

const gaussianCvar99 = gaussianMean + 2.33 * gaussianStddev;
const gaussianVsPareto = gaussianCvar99 > 0 ? Math.round((cvar99 / gaussianCvar99) * 1000) / 1000 : 0;

F7. Antifragility Score#

Scores each node on a [-1, +1] scale based on its IEC 62443 Security Level Target.

File: mc-engine.ts:791-795

Mathematical notation:

Scores: SL-T=0 -> -1.0 (fragile), SL-T=2 -> 0.0 (neutral), SL-T=4 -> +1.0 (antifragile).

F8. Barbell Score#

Measures defence budget concentration using coefficient of variation of cost distribution across layers. Range [0, 1].

File: mc-engine.ts:814-821

Mathematical notation:

const coV = lcMean > 0 ? lcStddev / lcMean : 0;
barbellScore = Math.round(Math.min(1, coV / 2) * 1000) / 1000;

F9. Layer Transition Conditional Probability Table#

Modifies edge weight when crossing between CDT layers using pre-calibrated transition probabilities.

File: mc-engine.ts:196-201

const LAYER_CPT: Record<string, number> = {
    "CS1->CS2": 0.40, "CS2->CS0": 0.25, "CS0->CS3": 0.70, "CS1->CS3": 0.60,
    "CS3->CS5": 0.85, "CS5->CS4": 0.15, "CS1->CS6": 0.30, "CS3->CS7": 0.50,
    "CS1->CS0": 0.35, "CS2->CS3": 0.55, "CS4->CS1": 0.45, "CS6->CS1": 0.40,
    "CS7->CS5": 0.35, "CS0->CS2": 0.50, "CS3->CS4": 0.20, "CS5->CS7": 0.30,
};
// Edge weight adjusted: w_adjusted = w_raw * CPT[src_layer -> tgt_layer]
// Default 0.3 for unmapped transitions

3. Edge Weight Composition#

F10. 14-Dimension Edge Weight#

The core weight function combines a base weight per relationship type with 14 modifier dimensions from node properties and temporal signals.

File: mc-weights.ts:137-186

Mathematical notation:

Clamped to [0.01, 1.0].

export function computeEdgeWeight(relType: string, sourceProps: NodeProps, targetProps: NodeProps): number {
    let w = BASE_WEIGHTS[relType] ?? 0.25;
    if (targetProps.epss_score != null && targetProps.epss_score > 0)
        w *= (0.3 + 0.7 * targetProps.epss_score);
    if (targetProps.cvss_v3_score != null && targetProps.cvss_v3_score > 0)
        w *= (0.5 + 0.5 * (targetProps.cvss_v3_score / 10));
    if (targetProps.kev_listed) w *= 1.5;
    if (targetProps.sl_target != null && targetProps.sl_target > 0)
        w *= (1 - targetProps.sl_target * 0.18);
    if (sourceProps.tacam_score != null && sourceProps.tacam_score > 0)
        w *= (0.5 + 0.5 * sourceProps.tacam_score);
    if (sourceProps.erika_activation != null)
        w *= Math.max(0.2, sourceProps.erika_activation);
    // ... temporal signals (F8-F11) ...
    return Math.min(1.0, Math.max(0.01, w));
}

F11. EPSS Velocity Boost Map#

CVEs with rising EPSS scores (delta_30d > 0.05) get boosted during walk traversal.

File: mc-weights.ts:239-248

Range: [1.15, 1.5]. Source: seldon.epss_trajectory (555K rows), cached 10 minutes.

F12. TACAM Actor Recency Modifier#

Maps campaign recency score to a weight modifier for actor-sourced edges.

File: mc-weights.ts:296-304

Range: recency 0.0 (dormant) -> 0.3x, recency 1.0 (active) -> 1.5x.

F13. Spectral Vulnerability Boost#

Maps eigenvector centrality rank to a walk probability boost for critical pivot nodes.

File: mc-engine.ts:174-183

Range: eigen_rank=0 -> 1.8x, eigen_rank=0.1 -> 1.1x. Source: seldon.spectral_analysis, top 50 nodes with eigen_rank < 0.1.

4. Hawkes Self-Exciting Process#

F14. Hawkes Intensity Function#

The conditional intensity at time t given event history {t_i}.

File: mc-hawkes.ts:74-88

Mathematical notation:

Where:

• mu = background intensity (default 0.05, ~1 event per 20 hours)
• alpha = excitation amplitude (default 0.8)
• beta = decay rate (default 0.3, half-life ~2.3 hours)

function hawkesIntensity(t: number, eventTimes: number[], mu: number, alpha: number, beta: number): number {
    let lambda = mu;
    for (const ti of eventTimes) {
        if (ti < t) {
            lambda += alpha * Math.exp(-beta * (t - ti));
        }
    }
    return lambda;
}

F15. Basic Reproduction Number (R0)#

The expected number of secondary infections per primary infection. If R0 > 1, the cascade grows.

File: mc-hawkes.ts:294

Default: 0.8 / 0.3 = 2.67 (supercritical ; cascades grow).

F16. Generalized Pareto Distribution (GPD) Severity#

Samples loss severity from a heavy-tailed distribution.

File: mc-hawkes.ts:96-102

Mathematical notation:

For shape parameter xi > 0:

For xi = 0 (exponential limit):

function gpdSample(xi: number, sigma: number): number {
    const u = Math.random();
    if (Math.abs(xi) < 1e-8) {
        return -sigma * Math.log(1 - u); // exponential limit
    }
    return (sigma / xi) * (Math.pow(1 - u, -xi) - 1);
}

Default: xi = 0.6 (heavy tail), sigma = 2.0 ($M scale).

F17. Ogata Thinning Algorithm#

Efficient event simulation via acceptance-rejection on the Hawkes intensity upper bound.

File: mc-hawkes.ts:110-136

1. Compute upper bound: lambda_bar = lambda(t) + 0.01
2. Generate candidate inter-event time: w = -ln(U) / lambda_bar
3. Accept with probability: lambda(t + w) / lambda_bar

F18. Hawkes-SIR Cost with Intensity Boost#

GPD-sampled cost is amplified by the current Hawkes intensity relative to baseline.

File: mc-hawkes.ts:251-252

F19. SIR Recovery Model#

Infected nodes recover with probability following an exponential CDF with rate 1/72 (72-hour MTTR).

File: mc-hawkes.ts:197-205

5. Eight-Stage Structural Causal Model#

Pearl's do-calculus implemented as stage-specific causal mechanisms. Each stage receives upstream variables and propagates downstream with exogenous noise U ~ Uniform(0,1).

F20. CS0 ; Physical Process (Bernoulli Flow)#

File: mc-scm.ts:50-69

F21. CS1 ; Cyber Detection (EPSS-Calibrated)#

File: mc-scm.ts:71-93

F22. CS2 ; OT/ICS Isolation (IEC 62443 Zone Model)#

File: mc-scm.ts:95-116

F23. CS3 ; Organizational Impact (Business Continuity)#

File: mc-scm.ts:118-142

F24. CS4 ; Geopolitical Propagation (Leontief Input-Output)#

File: mc-scm.ts:144-174

Where g_risk is the mean g_risk_multiplier from public.governance_risk.

F25. CS5 ; Economic Quantification (ALE Model)#

File: mc-scm.ts:176-207

All values in $M. Premium multiplier: 1.3 + P_cascade * 0.7. Insurance coverage: 60 percent typical.

F26. CS6 ; Psychographic Shift (ERIKA Quantum State)#

File: mc-scm.ts:209-230

F27. CS7 ; Temporal Forecast (Seldon Prediction)#

File: mc-scm.ts:232-254

6. ALE Insurance Engine#

F28. Poisson-Pareto ALE (Core Formula)#

Annual Loss Expectancy via compound Poisson-Pareto Monte Carlo.

File: ale-engine.ts:894-915

Mathematical notation:

for (let i = 0; i < N_SIMS; i++) {
    // Poisson sampling (inverse CDF)
    let incidents = 0;
    let p = Math.exp(-adjustedFrequency);
    let cumP = p;
    const u = rng();
    while (u > cumP) {
        incidents++;
        p *= adjustedFrequency / incidents;
        cumP += p;
    }
    let totalLoss = 0;
    for (let j = 0; j < incidents; j++) {
        let loss = samplePareto(rng, sectorAlpha, sectorXmin);
        if (coverage_limit !== undefined && loss > coverage_limit) loss = coverage_limit;
        totalLoss += loss;
    }
    annualLosses[i] = totalLoss;
}

Default: 50,000 simulations. 17 sector-specific frequency calibrations.

F29. Frequency Adjustment (SL-T + Posture)#

File: ale-engine.ts:791-794

F30. Percentile-Based Tail Shape Estimator (Pareto Alpha)#

Robust alternative to the Hill estimator for mixed-magnitude datasets.

File: ale-engine.ts:830-833

Clamped to [1.5, 3.0] to guarantee finite variance.

F31. Reporting Bias Correction (xmin)#

Corrects for the fact that the threat_incidents table contains only publicly reported major incidents.

File: ale-engine.ts:854-857

Where R is annual revenue and 0.25 is the reporting bias factor (Romanosky 2016, Advisen).

Revenue scaling:

F32. Kolmogorov-Smirnov Goodness-of-Fit Test#

Tests whether the fitted Pareto distribution is consistent with the empirical tail, meaning the observed loss sample itself rather than a known ground truth distribution.

The word carries its statistical sense on this line and no wider one. $F_{emp}$ is the empirical distribution function of the sample drawn from the threat_incidents table, and the statistic compares the fitted curve against that sample alone. A small $D_n$ says the sample does not refute the Pareto fit at the chosen significance level, independently of the reporting-bias correction the preceding formula applies to $x_{min}$. The formula is stated here as part of the model specification, not as a result.

File: ale-engine.ts:225-268

Where:

P-value approximation (Marsaglia et al. 2003, simplified with Stephens correction):

Good fit if D_n < 1.36 / sqrt(n) at alpha=0.05.

F33. Gordon-Loeb Optimal Security Investment#

Maximum economically rational security spend.

File: ale-engine.ts:978

Based on Gordon & Loeb (2002): optimal investment never exceeds 1/e of expected loss.

F34. Premium Calculation#

Insurance premium with risk and expense loading.

File: ale-engine.ts:967-975

Where:

• 1.25 = 25 percent risk margin
• 1.15 = 15 percent admin/acquisition costs
• Floor = 0.03 percent of revenue (minimum market premium)

90 percent confidence interval: [P5 * 1.4375, P95 * 1.4375].

F35. Y5381 War/Terrorism Exclusion#

Lloyd's Market Bulletin Y5381 state-backed loss attribution.

File: ale-engine.ts:706-730

Where state_portion = cost-weighted fraction of incidents attributed to state-sponsored actors (from seldon.actor_eic attribution + public.threat_incidents.nation_state_suspected). Clamped to [0, 0.95].

War exclusion applicable if state_portion > 0.10.

F36. Loss Development Factors (LDF)#

Fraction of ultimate loss reported at each development year, by sector.

File: ale-engine.ts:283-299

Sources: Advisen cyber loss data, NetDiligence claims studies (2020-2025).

F37. xoshiro128** PRNG (ALE Engine)#

Seeded PRNG for the ALE engine (different from Mulberry32 in mc-engine).

File: ale-engine.ts:366-379

function createRng(seed: number): () => number {
    let s0 = seed | 0 || 1;
    let s1 = (seed * 2654435761) | 0 || 1;
    let s2 = (seed * 2246822519) | 0 || 1;
    let s3 = (seed * 3266489917) | 0 || 1;
    return () => {
        const result = (((s1 * 5) << 7 | (s1 * 5) >>> 25) * 9) >>> 0;
        const t = s1 << 9;
        s2 ^= s0; s3 ^= s1; s1 ^= s2; s0 ^= s3;
        s2 ^= t;
        s3 = s3 << 11 | s3 >>> 21;
        return result / 4294967296;
    };
}

7. ATQ Actor Threat Quotient#

8-component weighted scoring system computed entirely in Postgres stored procedures. Output range: 0-100 via sigmoid transformation.

F38. ATQ Weighted Sum and Sigmoid#

File: atq-migration.sql:598-612

The intercept -4.0 centres the sigmoid: max Z = 8.0 (all C=1.0), so Z in [-4.0, +4.0]. Median actors (~sum 4.0) map to ATQ ~50; top actors (~sum 6.0) map to ATQ ~88.

v_z := 1.8 * v_c1 + 1.4 * v_c2 + 1.2 * v_c3 + 1.0 * v_c4
     + 0.8 * v_c5 + 0.7 * v_c6 + 0.6 * v_c7 + 0.5 * v_c8
     - 4.0;
v_atq := (1.0 / (1.0 + EXP(-v_z))) * 100.0;

F38a. ATQ Confidence Interval#

File: atq-migration.sql:614-624

Clamped to [0, 100].

Component Formulas (C1-C8)#

C1: EIC Composite (CS6 Psychographic)#

File: atq-migration.sql:178-208

Where I = intent, C = capability, O = opportunity (from actor_eic), D = Dark Triad d_factor (from psychometric_profiles).

C2: TACAM Affinity (CS1-CS2 Cyber-Physical)#

File: atq-migration.sql:213-258

Where:

• technique_breadth = actor_techniques / max_techniques (across all actors)
• sector_reach = distinct_sectors(score > 0.3) / 17
• protocol_reach = distinct_protocols(score > 0.3) / 12
• kill_chain_completeness = distinct_tactics / 14

C3: Temporal Momentum (CS7 Temporal)#

File: atq-migration.sql:283-337

Where:

• $R_{campaign} = \exp\left(-\frac{\ln 2 \times \Delta t}{90 \text{ days}}\right)$ (half-life 90 days)
• $V_{EPSS} = \text{clamp}\left(\frac{\text{mean\_epss\_delta\_30d}}{0.05}, 0, 1\right)$
• $V_{technique} = \min\left(\frac{\text{techniques\_added\_last\_year}}{10}, 1\right)$
• $R_{incident} = \exp\left(-\frac{\ln 2 \times \Delta t}{180 \text{ days}}\right)$ (half-life 180 days)

C4: Incident Evidence (CS5 Economic)#

File: atq-migration.sql:342-381

Where:

• $F = \min\left(\frac{\text{count}}{20}, 1\right)$ (incidents in last 2 years)
• $S = \min(1, \ \text{pct}_{crit} \times 1.0 + \text{pct}_{high} \times 0.7 + \text{pct}_{med} \times 0.3)$
• $K = \min\left(\frac{\log_{10}(\text{total\_cost\_m} + 1)}{4}, 1\right)$

C5: Exploit Economics Index (CS5/CS7)#

File: atq-migration.sql:386-427

With direct EEI data: $C_5 = \min(1, \ 0.60 \times EEI + 0.40 \times EPSS_{avg})$

With EPSS proxy: $C_5 = \min(1, \ EPSS_{avg} \times 1.5 + \frac{KEV_{count}}{CWE_{count}} \times 0.2)$

C6: Discourse Dynamics (CS6 Psychographic)#

File: atq-migration.sql:431-468

Where:

• D_base: Hysteric=0.80, Master=0.65, Analyst=0.60, University=0.40, else=0.50
• B_bifurcation = 0.15 if beta_tension < 0.15 or near_bifurcation = true, else 0
• M_stability = 0.10 * (1 - beta_tension)

C7: Geopolitical Pressure (CS4 Geopolitical)#

File: atq-migration.sql:472-511

Where:

• T_origin = max(geo_risk_score * attribution_weight) for attribution > 0.3
• T_target = avg(geo_risk_score * target_weight) for target > 0.1
• C_conflict = (conflict_exposure + state_sponsorship_likelihood) / 2 from geopolitical_field

C8#

Kramers Barrier Penetration (CS0-CS2 Physics)

File: atq-migration.sql:515-561

Where h = barrier height, h_max = global max, and f_adj is a discourse-adjusted factor:

Low barrier = high penetration = high score (inverse relationship).

8. Detection Model#

F39. IEC 62443 Zone-Aware Detection Probability#

Per-hop detection probability during Monte Carlo walks, accounting for zone transitions and SL-T gaps.

File: mc-engine.ts:598-601

Where:

• Z_cross = 1.3 if zone transition, 1.0 otherwise
• G_penalty = 0.05 if SL-T = 0 (unprotected zone, 95 percent detection reduction), 1.0 otherwise

const gapPenalty = slTarget === 0 ? 0.05 : 1.0;
const pDetected = Math.min(0.95,
    (slTarget / 4) * 0.5 * (isZoneTransition ? 1.3 : 1.0) * gapPenalty
);

9. Spectral Vulnerability#

See F13 above. The spectral boost map is loaded from seldon.spectral_analysis (refreshed every 15 minutes). Nodes with eigen_rank < 0.1 (top 10 percent eigenvector centrality) are the most critical graph pivot points.

10. EPSS Velocity#

See F11 above. EPSS velocity data is sourced from seldon.epss_trajectory (555,556 rows). CVEs with delta_30d > 0.05 are flagged. The boost is applied both during edge weight computation (F10, dimension 8) and as a per-walk overlay (F11, dimension 13). Cache TTL: 10 minutes.

11. Formula Index#

The table below indexes every formula in this reference by its file and line, so a reader auditing the codebase against the mathematics can locate either from the other.

12. References#

This document catalogues formulas already implemented in the Cyber Digital Twin codebase rather than external research, and every formula above names its own source file and line, mc-engine.ts, mc-weights.ts, mc-hawkes.ts, mc-scm.ts, ale-engine.ts and atq-migration.sql, in place of a citation. No published source is cited in this reference.