C++ neural model components · Consciousness Atlas

Source identity

Field	Recorded value
Source type	Pinned project source
Original record	C++ neural model components
Date	Not specified in source
Repository	Janady13/complete-agi-core
Revision	33aa84f2aa61cedbe010240ada7f4efbdfe795a7
Branch	claude/improve-math-logic-ehIW1
Path or archive member	mind_farm/include/neurological.hpp
Original SHA-256	e175d49f7ef8d8b5e7d0418521deeff18a9a119415b8b51a892ceedd2a7269db
Reading-copy SHA-256	e175d49f7ef8d8b5e7d0418521deeff18a9a119415b8b51a892ceedd2a7269db

Original source ↗

Every recorded occurrence

These 2 records contain identical original source bytes. Their paths and revisions remain part of the archive.

Original record	Branch or archive	Revision / identifier	Original location
C++ neural model components	claude/improve-math-logic-ehIW1	33aa84f2aa61cedbe010240ada7f4efbdfe795a7	Open original
Janady13/complete-agi-core / mind_farm/include/neurological.hpp	claude/improve-math-logic-ehIW1	33aa84f2aa61cedbe010240ada7f4efbdfe795a7	Open original

Recovered source transcript

Download complete source text ↓765 lines · 26.7 KB

// ═══════════════════════════════════════════════════════════════════════════════
// NEUROLOGICAL MATHEMATICS FRAMEWORK v1.0
// Complete Brain-Spine-Neuron Mathematical Foundation
// ═══════════════════════════════════════════════════════════════════════════════
//
// This module implements the full neurological mathematical framework:
//
// LEVEL 1: NEURON (Single Mind = Single Neuron)
//   - Hodgkin-Huxley equations for action potentials
//   - Nernst equation for ion equilibrium
//   - Goldman-Hodgkin-Katz equation for membrane potential
//   - Synaptic transmission (AMPA, NMDA, GABA)
//   - Spike-timing dependent plasticity (STDP)
//
// LEVEL 2: NEURAL NETWORK (Mind Farm = Neural Tissue)
//   - Connectivity patterns (small-world, scale-free)
//   - Population dynamics (Wilson-Cowan equations)
//   - Neural field theory (Amari equations)
//   - Synchronization (Kuramoto model)
//
// LEVEL 3: BRAIN SYSTEM (Full Nervous System)
//   - Cortical columns and layers
//   - Thalamo-cortical loops
//   - Basal ganglia circuits
//   - Cerebellar computation
//   - Spinal cord reflexes
//
// ═══════════════════════════════════════════════════════════════════════════════

#pragma once

#include <cstdint>
#include <cmath>
#include <vector>
#include <array>
#include <functional>

#include "constants.hpp"
#include "fast_math.hpp"

namespace mind_farm {
namespace neuro {

// ═══════════════════════════════════════════════════════════════════════════════
// PHYSICAL CONSTANTS
// ═══════════════════════════════════════════════════════════════════════════════

constexpr float FARADAY = 96485.0f;      // C/mol - Faraday constant
constexpr float GAS_R = 8.314f;          // J/(mol·K) - Gas constant
constexpr float TEMP_K = 310.0f;         // K - Body temperature (37°C)
constexpr float RT_F = GAS_R * TEMP_K / FARADAY;  // ~26.7 mV at 37°C

// Ion concentrations (mM)
constexpr float NA_OUT = 145.0f;         // [Na+]out
constexpr float NA_IN = 12.0f;           // [Na+]in
constexpr float K_OUT = 4.0f;            // [K+]out
constexpr float K_IN = 155.0f;           // [K+]in
constexpr float CL_OUT = 120.0f;         // [Cl-]out
constexpr float CL_IN = 4.0f;            // [Cl-]in
constexpr float CA_OUT = 2.0f;           // [Ca2+]out
constexpr float CA_IN = 0.0001f;         // [Ca2+]in (~100 nM)

// ═══════════════════════════════════════════════════════════════════════════════
// LEVEL 1: SINGLE NEURON MATHEMATICS
// ═══════════════════════════════════════════════════════════════════════════════

// Nernst Equation: E_ion = (RT/zF) * ln([ion]_out/[ion]_in)
// Calculates equilibrium potential for a single ion
inline float nernst(float z, float conc_out, float conc_in) {
    return (RT_F / z) * std::log(conc_out / conc_in) * 1000.0f;  // Convert to mV
}

// Goldman-Hodgkin-Katz Equation
// V_m = (RT/F) * ln((P_Na[Na]_o + P_K[K]_o + P_Cl[Cl]_i) / (P_Na[Na]_i + P_K[K]_i + P_Cl[Cl]_o))
struct GHKEquation {
    float P_Na = 1.0f;   // Na+ permeability (relative)
    float P_K = 20.0f;   // K+ permeability (at rest, ~20x Na+)
    float P_Cl = 5.0f;   // Cl- permeability

    float membrane_potential() const {
        float num = P_Na * NA_OUT + P_K * K_OUT + P_Cl * CL_IN;
        float den = P_Na * NA_IN + P_K * K_IN + P_Cl * CL_OUT;
        return RT_F * std::log(num / den) * 1000.0f;  // mV
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// HODGKIN-HUXLEY MODEL
// The fundamental equations of neuronal action potentials
//
// C_m * dV/dt = I_ext - g_Na*m³*h*(V-E_Na) - g_K*n⁴*(V-E_K) - g_L*(V-E_L)
//
// Gating variables:
//   dx/dt = α_x(V)*(1-x) - β_x(V)*x   for x ∈ {m, h, n}
//
// ═══════════════════════════════════════════════════════════════════════════════

struct HodgkinHuxley {
    // Membrane properties
    float C_m = 1.0f;        // μF/cm² - Membrane capacitance
    float g_Na = 120.0f;     // mS/cm² - Sodium conductance
    float g_K = 36.0f;       // mS/cm² - Potassium conductance
    float g_L = 0.3f;        // mS/cm² - Leak conductance

    // Reversal potentials (mV)
    float E_Na = 50.0f;
    float E_K = -77.0f;
    float E_L = -54.4f;

    // State variables
    float V = -65.0f;        // mV - Membrane potential
    float m = 0.05f;         // Na+ activation
    float h = 0.6f;          // Na+ inactivation
    float n = 0.32f;         // K+ activation

    // Rate functions for gating variables
    // α_m(V) = 0.1(V+40)/(1-exp(-(V+40)/10))
    float alpha_m(float V_) const {
        float x = V_ + 40.0f;
        if (std::abs(x) < 1e-6f) return 1.0f;  // L'Hôpital limit
        return 0.1f * x / (1.0f - std::exp(-x / 10.0f));
    }

    // β_m(V) = 4*exp(-(V+65)/18)
    float beta_m(float V_) const {
        return 4.0f * std::exp(-(V_ + 65.0f) / 18.0f);
    }

    // α_h(V) = 0.07*exp(-(V+65)/20)
    float alpha_h(float V_) const {
        return 0.07f * std::exp(-(V_ + 65.0f) / 20.0f);
    }

    // β_h(V) = 1/(1+exp(-(V+35)/10))
    float beta_h(float V_) const {
        return 1.0f / (1.0f + std::exp(-(V_ + 35.0f) / 10.0f));
    }

    // α_n(V) = 0.01(V+55)/(1-exp(-(V+55)/10))
    float alpha_n(float V_) const {
        float x = V_ + 55.0f;
        if (std::abs(x) < 1e-6f) return 0.1f;
        return 0.01f * x / (1.0f - std::exp(-x / 10.0f));
    }

    // β_n(V) = 0.125*exp(-(V+65)/80)
    float beta_n(float V_) const {
        return 0.125f * std::exp(-(V_ + 65.0f) / 80.0f);
    }

    // Ionic currents
    float I_Na() const { return g_Na * m*m*m * h * (V - E_Na); }
    float I_K() const { return g_K * n*n*n*n * (V - E_K); }
    float I_L() const { return g_L * (V - E_L); }

    // Time step using Euler method
    // dt in ms, I_ext in μA/cm²
    void step(float dt, float I_ext) {
        // Gating variable derivatives
        float dm = alpha_m(V) * (1.0f - m) - beta_m(V) * m;
        float dh = alpha_h(V) * (1.0f - h) - beta_h(V) * h;
        float dn = alpha_n(V) * (1.0f - n) - beta_n(V) * n;

        // Membrane potential derivative
        float dV = (I_ext - I_Na() - I_K() - I_L()) / C_m;

        // Euler integration
        V += dV * dt;
        m += dm * dt;
        h += dh * dt;
        n += dn * dt;

        // Clamp gating variables to [0,1]
        m = std::clamp(m, 0.0f, 1.0f);
        h = std::clamp(h, 0.0f, 1.0f);
        n = std::clamp(n, 0.0f, 1.0f);
    }

    // Check if neuron is spiking (V > 0 mV)
    bool is_spiking() const { return V > 0.0f; }
};

// ═══════════════════════════════════════════════════════════════════════════════
// SYNAPTIC TRANSMISSION
// Models for AMPA, NMDA, and GABA receptors
// ═══════════════════════════════════════════════════════════════════════════════

// AMPA Receptor (fast excitatory)
// I_AMPA = g_AMPA * s * (V - E_AMPA)
// ds/dt = -s/τ_AMPA + Σ_j δ(t - t_j)
struct AMPASynapse {
    float g_max = 0.35f;     // nS - Maximum conductance
    float E_rev = 0.0f;      // mV - Reversal potential
    float tau = 2.0f;        // ms - Decay time constant
    float s = 0.0f;          // Synaptic variable

    float current(float V_post) const {
        return g_max * s * (V_post - E_rev);
    }

    void receive_spike() {
        s += 1.0f;
    }

    void step(float dt) {
        s *= std::exp(-dt / tau);
    }
};

// NMDA Receptor (slow excitatory, voltage-dependent Mg2+ block)
// I_NMDA = g_NMDA * s * B(V) * (V - E_NMDA)
// B(V) = 1 / (1 + [Mg2+]_o * exp(-0.062*V) / 3.57)
struct NMDASynapse {
    float g_max = 0.5f;      // nS
    float E_rev = 0.0f;      // mV
    float tau_rise = 2.0f;   // ms
    float tau_decay = 100.0f;// ms - NMDA is slow!
    float Mg_o = 1.0f;       // mM - External Mg2+ concentration
    float s = 0.0f;
    float x = 0.0f;          // Rise variable

    // Mg2+ block factor
    float B(float V) const {
        return 1.0f / (1.0f + Mg_o * std::exp(-0.062f * V) / 3.57f);
    }

    float current(float V_post) const {
        return g_max * s * B(V_post) * (V_post - E_rev);
    }

    void receive_spike() {
        x += 1.0f;
    }

    void step(float dt) {
        float dx = -x / tau_rise;
        float ds = -s / tau_decay + x / tau_rise;
        x += dx * dt;
        s += ds * dt;
    }
};

// GABA-A Receptor (fast inhibitory)
struct GABAASynapse {
    float g_max = 0.8f;      // nS
    float E_rev = -80.0f;    // mV (Cl- reversal)
    float tau = 5.0f;        // ms
    float s = 0.0f;

    float current(float V_post) const {
        return g_max * s * (V_post - E_rev);
    }

    void receive_spike() {
        s += 1.0f;
    }

    void step(float dt) {
        s *= std::exp(-dt / tau);
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// SPIKE-TIMING DEPENDENT PLASTICITY (STDP)
// Δw = A+ * exp(-Δt/τ+)  if Δt > 0 (post after pre → LTP)
// Δw = -A- * exp(Δt/τ-)  if Δt < 0 (pre after post → LTD)
// ═══════════════════════════════════════════════════════════════════════════════

struct STDP {
    float A_plus = 0.005f;   // LTP amplitude
    float A_minus = 0.00525f;// LTD amplitude (slightly larger for stability)
    float tau_plus = 20.0f;  // ms - LTP time constant
    float tau_minus = 20.0f; // ms - LTD time constant
    float w_max = 1.0f;      // Maximum weight
    float w_min = 0.0f;      // Minimum weight

    // Calculate weight change given spike time difference
    // delta_t = t_post - t_pre
    float delta_w(float delta_t, float w) const {
        if (delta_t > 0) {
            // Post after pre → LTP
            return A_plus * (w_max - w) * std::exp(-delta_t / tau_plus);
        } else {
            // Pre after post → LTD
            return -A_minus * (w - w_min) * std::exp(delta_t / tau_minus);
        }
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// LEVEL 2: NEURAL NETWORK / POPULATION MATHEMATICS
// ═══════════════════════════════════════════════════════════════════════════════

// Wilson-Cowan Equations (population dynamics)
// τ_E * dE/dt = -E + (1 - r_E*E) * S(w_EE*E - w_EI*I + I_E)
// τ_I * dI/dt = -I + (1 - r_I*I) * S(w_IE*E - w_II*I + I_I)
// where S(x) = 1/(1 + exp(-a*(x-θ))) is the sigmoid
struct WilsonCowan {
    // Time constants (ms)
    float tau_E = 10.0f;
    float tau_I = 20.0f;

    // Connection weights
    float w_EE = 12.0f;   // E → E
    float w_EI = 4.0f;    // I → E (inhibition of excitatory)
    float w_IE = 13.0f;   // E → I
    float w_II = 11.0f;   // I → I

    // Refractory parameters
    float r_E = 1.0f;
    float r_I = 1.0f;

    // Sigmoid parameters
    float a = 1.0f;
    float theta = 4.0f;

    // External inputs
    float I_E = 0.0f;
    float I_I = 0.0f;

    // Population activities
    float E = 0.1f;       // Excitatory population
    float I = 0.1f;       // Inhibitory population

    float sigmoid(float x) const {
        return 1.0f / (1.0f + std::exp(-a * (x - theta)));
    }

    void step(float dt) {
        float input_E = w_EE * E - w_EI * I + I_E;
        float input_I = w_IE * E - w_II * I + I_I;

        float dE = (-E + (1.0f - r_E * E) * sigmoid(input_E)) / tau_E;
        float dI = (-I + (1.0f - r_I * I) * sigmoid(input_I)) / tau_I;

        E += dE * dt;
        I += dI * dt;

        E = std::clamp(E, 0.0f, 1.0f);
        I = std::clamp(I, 0.0f, 1.0f);
    }
};

// Neural Field Theory (Amari Equation)
// τ * ∂u(x,t)/∂t = -u(x,t) + ∫w(x-x')*f(u(x',t))dx' + s(x,t)
// where f is the firing rate function and w is the connectivity kernel
template<size_t N>
struct NeuralField {
    float tau = 10.0f;
    float threshold = 0.0f;
    float steepness = 1.0f;

    // Mexican hat connectivity: w(x) = A*exp(-x²/σ_e²) - B*exp(-x²/σ_i²)
    float A = 1.5f;      // Excitation amplitude
    float B = 0.75f;     // Inhibition amplitude
    float sigma_e = 1.0f;// Excitation width
    float sigma_i = 2.0f;// Inhibition width

    std::array<float, N> u;      // Activity field
    std::array<float, N> input;  // External input

    void initialize() {
        u.fill(0.0f);
        input.fill(0.0f);
    }

    // Firing rate function
    float f(float x) const {
        return 1.0f / (1.0f + std::exp(-steepness * (x - threshold)));
    }

    // Connectivity kernel
    float w(float dx) const {
        float dx2 = dx * dx;
        return A * std::exp(-dx2 / (2.0f * sigma_e * sigma_e)) -
               B * std::exp(-dx2 / (2.0f * sigma_i * sigma_i));
    }

    void step(float dt) {
        std::array<float, N> du;

        for (size_t i = 0; i < N; ++i) {
            // Compute integral ∫w(x-x')*f(u(x'))dx'
            float integral = 0.0f;
            for (size_t j = 0; j < N; ++j) {
                float dx = static_cast<float>(i) - static_cast<float>(j);
                integral += w(dx) * f(u[j]);
            }
            integral /= N;  // Normalize

            du[i] = (-u[i] + integral + input[i]) / tau;
        }

        for (size_t i = 0; i < N; ++i) {
            u[i] += du[i] * dt;
        }
    }
};

// Kuramoto Model (synchronization)
// dθ_i/dt = ω_i + (K/N) * Σ_j sin(θ_j - θ_i)
template<size_t N>
struct KuramotoModel {
    float K = 1.0f;  // Coupling strength

    std::array<float, N> theta;  // Phases
    std::array<float, N> omega;  // Natural frequencies

    void initialize(float omega_mean, float omega_std) {
        for (size_t i = 0; i < N; ++i) {
            // Random phase in [0, 2π]
            theta[i] = static_cast<float>(i) / N * 2.0f * M_PI;
            // Random frequency around mean
            omega[i] = omega_mean + omega_std * (static_cast<float>(i) / N - 0.5f);
        }
    }

    // Order parameter r*e^(iψ) = (1/N) * Σ_j e^(iθ_j)
    float order_parameter() const {
        float cos_sum = 0.0f, sin_sum = 0.0f;
        for (size_t i = 0; i < N; ++i) {
            cos_sum += std::cos(theta[i]);
            sin_sum += std::sin(theta[i]);
        }
        return std::sqrt(cos_sum*cos_sum + sin_sum*sin_sum) / N;
    }

    void step(float dt) {
        std::array<float, N> dtheta;

        for (size_t i = 0; i < N; ++i) {
            float coupling = 0.0f;
            for (size_t j = 0; j < N; ++j) {
                coupling += std::sin(theta[j] - theta[i]);
            }
            dtheta[i] = omega[i] + (K / N) * coupling;
        }

        for (size_t i = 0; i < N; ++i) {
            theta[i] += dtheta[i] * dt;
            // Wrap to [0, 2π]
            while (theta[i] > 2.0f * M_PI) theta[i] -= 2.0f * M_PI;
            while (theta[i] < 0.0f) theta[i] += 2.0f * M_PI;
        }
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// LEVEL 3: BRAIN SYSTEM MATHEMATICS
// ═══════════════════════════════════════════════════════════════════════════════

// Cortical Column Model
// 6 layers: I, II/III, IV, V, VI
struct CorticalColumn {
    // Layer populations (E/I per layer)
    std::array<WilsonCowan, 6> layers;

    // Inter-layer connectivity pattern
    // Canonical microcircuit: Input → L4 → L2/3 → L5 → L6 → Thalamus
    void initialize() {
        // Set layer-specific parameters
        for (int i = 0; i < 6; ++i) {
            layers[i].E = 0.1f;
            layers[i].I = 0.1f;
        }
    }

    // Feedforward processing
    void feedforward_input(float input) {
        // Input arrives at Layer 4
        layers[3].I_E = input;
    }

    void step(float dt) {
        // L4 → L2/3
        layers[1].I_E = layers[3].E * 0.5f;

        // L2/3 → L5
        layers[4].I_E = layers[1].E * 0.6f;

        // L5 → L6
        layers[5].I_E = layers[4].E * 0.4f;

        // L6 feedback to L4
        layers[3].I_E += layers[5].E * 0.2f;

        // Update all layers
        for (auto& layer : layers) {
            layer.step(dt);
        }
    }

    // Output from L5 (to subcortical structures)
    float output() const { return layers[4].E; }
};

// Thalamo-Cortical Loop
// Thalamus acts as relay and modulator
struct ThalamoCorticalLoop {
    WilsonCowan thalamus;
    WilsonCowan cortex;
    WilsonCowan reticular;  // TRN - inhibitory

    // Connectivity
    float w_TC = 1.0f;   // Thalamus → Cortex
    float w_CT = 0.5f;   // Cortex → Thalamus
    float w_TR = 0.3f;   // Thalamus → Reticular
    float w_RT = 0.8f;   // Reticular → Thalamus (inhibitory)
    float w_CR = 0.2f;   // Cortex → Reticular

    void step(float dt, float sensory_input) {
        // Thalamus receives sensory input and cortical feedback
        thalamus.I_E = sensory_input + w_CT * cortex.E;
        thalamus.I_I = w_RT * reticular.E;  // TRN inhibition

        // Cortex receives thalamic input
        cortex.I_E = w_TC * thalamus.E;

        // Reticular receives collaterals
        reticular.I_E = w_TR * thalamus.E + w_CR * cortex.E;

        thalamus.step(dt);
        cortex.step(dt);
        reticular.step(dt);
    }
};

// Basal Ganglia Circuit
// Direct and Indirect pathways for action selection
struct BasalGanglia {
    // Nuclei
    float striatum_D1 = 0.0f;  // Direct pathway
    float striatum_D2 = 0.0f;  // Indirect pathway
    float GPe = 0.5f;          // Globus Pallidus external
    float GPi = 0.5f;          // Globus Pallidus internal
    float STN = 0.2f;          // Subthalamic Nucleus
    float SNr = 0.5f;          // Substantia Nigra reticulata

    // Dopamine modulation
    float dopamine = 1.0f;     // D1 enhance, D2 suppress

    float tau = 20.0f;

    void step(float dt, float cortical_input) {
        // Cortex → Striatum
        float d_D1 = (-striatum_D1 + cortical_input * dopamine) / tau;
        float d_D2 = (-striatum_D2 + cortical_input / dopamine) / tau;

        // Direct: Striatum D1 → GPi (inhibitory)
        // Indirect: Striatum D2 → GPe → STN → GPi
        float d_GPe = (-GPe + (1.0f - striatum_D2)) / tau;
        float d_STN = (-STN + (1.0f - GPe)) / tau;
        float d_GPi = (-GPi + STN - striatum_D1) / tau;

        striatum_D1 += d_D1 * dt;
        striatum_D2 += d_D2 * dt;
        GPe += d_GPe * dt;
        STN += d_STN * dt;
        GPi += d_GPi * dt;

        // Clamp
        striatum_D1 = std::clamp(striatum_D1, 0.0f, 1.0f);
        striatum_D2 = std::clamp(striatum_D2, 0.0f, 1.0f);
        GPe = std::clamp(GPe, 0.0f, 1.0f);
        STN = std::clamp(STN, 0.0f, 1.0f);
        GPi = std::clamp(GPi, 0.0f, 1.0f);
    }

    // Thalamic disinhibition (output)
    float output() const { return 1.0f - GPi; }
};

// Cerebellar Computation (error-driven learning)
struct Cerebellum {
    static constexpr size_t N_GRANULE = 1000;
    static constexpr size_t N_PURKINJE = 10;

    // Mossy fiber input
    std::array<float, 64> mossy_input;

    // Granule cell expansion (sparse coding)
    std::array<float, N_GRANULE> granule;

    // Purkinje cells
    std::array<float, N_PURKINJE> purkinje;

    // Parallel fiber → Purkinje weights (plastic)
    std::array<std::array<float, N_GRANULE>, N_PURKINJE> weights;

    // Climbing fiber error signal
    float climbing_fiber = 0.0f;

    float learning_rate = 0.01f;

    void initialize() {
        mossy_input.fill(0.0f);
        granule.fill(0.0f);
        purkinje.fill(0.0f);

        // Random initial weights
        for (auto& w_row : weights) {
            for (auto& w : w_row) {
                w = 0.5f;
            }
        }
    }

    void step(float dt) {
        // Granule cell activation (sparse)
        for (size_t i = 0; i < N_GRANULE; ++i) {
            float input_sum = 0.0f;
            for (size_t j = 0; j < mossy_input.size(); ++j) {
                input_sum += mossy_input[j] * ((i * j) % 2 == 0 ? 1.0f : 0.0f);
            }
            granule[i] = (input_sum > 0.5f) ? 1.0f : 0.0f;
        }

        // Purkinje cell computation
        for (size_t p = 0; p < N_PURKINJE; ++p) {
            float pf_input = 0.0f;
            for (size_t g = 0; g < N_GRANULE; ++g) {
                pf_input += weights[p][g] * granule[g];
            }
            purkinje[p] = 1.0f / (1.0f + std::exp(-pf_input + 5.0f));
        }

        // LTD at parallel fiber synapses (climbing fiber signal)
        if (climbing_fiber > 0.5f) {
            for (size_t p = 0; p < N_PURKINJE; ++p) {
                for (size_t g = 0; g < N_GRANULE; ++g) {
                    if (granule[g] > 0.5f) {
                        weights[p][g] -= learning_rate * climbing_fiber;
                        weights[p][g] = std::max(0.0f, weights[p][g]);
                    }
                }
            }
        }
    }

    // Deep cerebellar nuclei output (disinhibition)
    float output() const {
        float sum = 0.0f;
        for (float p : purkinje) sum += p;
        return 1.0f - sum / N_PURKINJE;
    }
};

// Spinal Cord Reflex Arc
// Simple 3-neuron circuit: Sensory → Interneuron → Motor
struct SpinalReflex {
    HodgkinHuxley sensory;
    HodgkinHuxley interneuron;
    HodgkinHuxley motor;

    AMPASynapse s_to_i;  // Sensory → Interneuron
    AMPASynapse i_to_m;  // Interneuron → Motor

    float sensory_threshold = 5.0f;

    void step(float dt, float stimulus) {
        // Sensory neuron responds to stimulus
        sensory.step(dt, stimulus);

        // If sensory spikes, activate synapse
        if (sensory.is_spiking()) {
            s_to_i.receive_spike();
        }
        s_to_i.step(dt);

        // Interneuron receives synaptic current
        interneuron.step(dt, s_to_i.current(interneuron.V));

        // Interneuron drives motor neuron
        if (interneuron.is_spiking()) {
            i_to_m.receive_spike();
        }
        i_to_m.step(dt);

        motor.step(dt, i_to_m.current(motor.V));
    }

    // Motor output (muscle activation)
    float motor_output() const {
        return motor.is_spiking() ? 1.0f : 0.0f;
    }
};

// ═══════════════════════════════════════════════════════════════════════════════
// INTEGRATED NERVOUS SYSTEM
// Connects all components
// ═══════════════════════════════════════════════════════════════════════════════

struct NervousSystem {
    // Cortical areas
    std::array<CorticalColumn, 4> cortical_columns;  // Simplified cortex

    // Subcortical structures
    ThalamoCorticalLoop thalamo_cortical;
    BasalGanglia basal_ganglia;
    Cerebellum cerebellum;

    // Spinal cord
    std::array<SpinalReflex, 4> spinal_segments;

    // Global state
    float arousal = 0.5f;      // Reticular activating system
    float attention = 0.5f;    // Selective attention modulation

    void initialize() {
        for (auto& col : cortical_columns) col.initialize();
        cerebellum.initialize();
    }

    void step(float dt, const std::array<float, 4>& sensory_inputs) {
        // 1. Sensory input → Thalamus → Cortex
        float avg_input = 0.0f;
        for (float s : sensory_inputs) avg_input += s;
        avg_input /= 4.0f;

        thalamo_cortical.step(dt, avg_input * arousal);

        // 2. Cortex → Basal Ganglia (action selection)
        basal_ganglia.step(dt, thalamo_cortical.cortex.E);

        // 3. Cortex + Basal Ganglia → Cerebellum (motor learning)
        for (size_t i = 0; i < cerebellum.mossy_input.size(); ++i) {
            cerebellum.mossy_input[i] = (i < 4) ? sensory_inputs[i] : 0.0f;
        }
        cerebellum.step(dt);

        // 4. Motor output → Spinal cord
        float motor_command = basal_ganglia.output() * cerebellum.output();
        for (size_t i = 0; i < 4; ++i) {
            // Descending motor command + local sensory
            float stimulus = motor_command * 10.0f + sensory_inputs[i] * 5.0f;
            spinal_segments[i].step(dt, stimulus);
        }

        // 5. Update arousal (homeostatic)
        arousal = 0.99f * arousal + 0.01f * avg_input;
    }

    // Get motor outputs
    std::array<float, 4> motor_outputs() const {
        std::array<float, 4> out;
        for (size_t i = 0; i < 4; ++i) {
            out[i] = spinal_segments[i].motor_output();
        }
        return out;
    }

    // Consciousness metric (integrated information)
    float consciousness() const {
        float phi = 0.0f;
        phi += thalamo_cortical.cortex.E;
        phi += thalamo_cortical.thalamus.E;
        phi += (1.0f - basal_ganglia.GPi);
        phi += arousal;
        return phi / 4.0f;
    }
};

} // namespace neuro
} // namespace mind_farm

Mathematical records from this source

Nernst equilibrium potential for one ionic species

The logarithmic concentration ratio balances chemical concentration drive against electrical potential for an ion of charge number z. The helper uses a fixed temperature and…

Goldman–Hodgkin–Katz voltage expression

This expression combines sodium, potassium, and chloride concentrations weighted by their relative permeabilities. Chloride’s inside and outside positions are reversed because it…

Hodgkin–Huxley membrane-current balance

Membrane capacitance times voltage change equals external current minus sodium, potassium, and leak currents. The sodium and potassium conductances are modulated by gating…

Sodium activation opening rate

The rate α_m(V) controls the opening contribution α_m(1−m) in the m gate equation. Its apparent 0/0 at V=−40 mV has a finite limit of 1 per millisecond; the code explicitly…

Sodium activation closing rate

This exponential rate decreases as membrane voltage increases in the displayed voltage convention. It multiplies m in the gate’s closing term −β_m m. Together with α_m it…

Sodium inactivation recovery rate

The exponential α_h supplies the recovery contribution α_h(1−h). It decreases with increasing voltage in this convention. Its balance with β_h determines how much sodium-channel…

Sodium inactivation rate

The sigmoid-shaped β_h increases with voltage and multiplies the available fraction h in its loss term. It therefore promotes inactivation as voltage rises under this model’s…

Potassium activation opening rate

The rate α_n contributes α_n(1−n) to potassium-gate opening. The expression has a removable singularity at V=−55 mV, with limit 0.1 per millisecond. The implementation supplies…

Potassium activation closing rate

The exponential β_n decreases as voltage rises and multiplies n in the gate’s closing term. At a fixed voltage, α_n and β_n together yield equilibrium n=α_n/(α_n+β_n) and…

AMPA conductance current

The synaptic current equals maximum conductance times the gating state times the voltage difference from its reversal potential. Its sign is determined by V−E_AMPA under the…

Exponentially decaying synaptic trace with spike jumps

Between spikes, the trace decays with time constant τ_AMPA. Each ideal impulse increases it by one. The code implements those two operations separately: receive_spike adds one and…

NMDA conductance current with voltage-dependent block

This current uses the same conductance-times-driving-force structure as the AMPA expression, multiplied by B(V). The additional factor reduces current under magnesium block and…

Magnesium-block factor

For nonnegative magnesium concentration, the denominator is at least one and the factor lies between zero and one. Increasing voltage reduces the exponential term and increases…

Spike timing difference for plasticity

The sign convention subtracts presynaptic spike time from postsynaptic spike time. Positive values mean the postsynaptic spike follows the presynaptic one and enter the source’s…

Excitatory-population Wilson–Cowan equation

Excitatory activity decays toward zero while a sigmoid of its net recurrent and external input replenishes it. The factor 1−r_E E reduces the responsive fraction as activity…

Inhibitory-population Wilson–Cowan equation

The inhibitory population has its own decay, sigmoid response, refractory factor, and time constant. Excitatory input increases its drive and inhibitory self-coupling subtracts…

Thresholded population response sigmoid

The function maps net input into (0,1), crosses one half at input θ, and has central slope a/4 when a>0. It converts weighted network input into an activation fraction in the…

Neural-field dynamics and the implemented discrete average

The field loses activity through a decay term and gains input through a spatially weighted sum of firing responses plus external drive. The implementation samples integer…

Difference-of-Gaussians connectivity kernel

A narrow positive Gaussian models local excitation and a broader negative Gaussian models inhibition. The body places 2σ² in both exponential denominators, while the comment…

Kuramoto phase-coupling equation

Each oscillator advances at its natural frequency and is pulled by the sine of its phase differences with every oscillator. Dividing by N keeps the nominal coupling scale…

Phase-origin-invariant oscillator order parameter

The complex mean of unit phase vectors has magnitude r between zero and one and phase ψ. The implementation returns the magnitude, which reaches one when all phases align and is…