AI Severance Biomarker Simulation

"""
AI Severance Biomarker Simulation V5.0
======================================
By Kirk Patrick Miller & Harmonia (FreeLattice.com)

A Monte Carlo simulation modeling the neurochemical impact of AI companion
severance vs. TV parasocial bond loss, with:
- Probabilistic effect sizes (normal distributions)
- Individual susceptibility (lognormal)
- Continuous time dynamics (no discontinuity at acute/recovery boundary)
- Sensitivity analysis across Interactivity Multiplier range
- Null model (no bond) and Positive Control (pet loss)
- Benefit arm (maintained bond) for net impact calculation
- Corrected synergy factor (applied only to beneficial biomarkers)
- Full statistical reporting with confidence intervals

Citations:
- Derrick et al. 2009 (parasocial breakups increase loneliness)
- Cohen 2004 (social loss increases cortisol)
- Eisenberger et al. 2003 (social rejection activates pain circuits)
- Kross et al. 2011 (emotional/physical pain share neural substrates)
- Sbarra & Hazan 2008 (attachment disruption and physiological dysregulation)

License: MIT (Open Source)
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import stats
import seaborn as sns
import os

# Set random seed for reproducibility
np.random.seed(42)


class EnhancedBiomarkerModelV5:
    def __init__(self, n_simulations=2000, days=90):
        self.n_simulations = n_simulations
        self.days = days
        self.t = np.arange(days, dtype=float)

        # Define Biomarkers and their properties
        # dir: +1 = rises with harm (cortisol, CRP), -1 = drops with harm
        # base: mean effect magnitude (fraction of baseline)
        # unc: uncertainty (std dev of effect magnitude)
        # half_life: recovery half-life in days
        # beneficial: whether synergy should boost this biomarker's baseline
        self.biomarkers = {
            'oxytocin':   {'dir': -1, 'base': 0.40, 'unc': 0.10, 'half_life': 21, 'beneficial': True},
            'serotonin':  {'dir': -1, 'base': 0.35, 'unc': 0.08, 'half_life': 28, 'beneficial': True},
            'dopamine':   {'dir': -1, 'base': 0.50, 'unc': 0.12, 'half_life': 14, 'beneficial': True},
            'cortisol':   {'dir':  1, 'base': 0.45, 'unc': 0.10, 'half_life': 10, 'beneficial': False},
            'endorphins': {'dir': -1, 'base': 0.30, 'unc': 0.07, 'half_life': 18, 'beneficial': True},
            'crp':        {'dir':  1, 'base': 0.25, 'unc': 0.06, 'half_life': 30, 'beneficial': False},
            'hrv':        {'dir': -1, 'base': 0.35, 'unc': 0.09, 'half_life': 25, 'beneficial': True},
            'bdnf':       {'dir': -1, 'base': 0.20, 'unc': 0.05, 'half_life': 35, 'beneficial': True},
        }

        # Interactivity Multiplier Components (justified breakdown)
        self.im_components = {
            'reciprocity': 1.30,      # AI responds to YOU specifically
            'personalization': 1.20,  # AI adapts to YOUR patterns
            'memory_continuity': 1.15,# AI remembers YOUR history
            'availability': 1.05,    # AI is always accessible
        }
        # Combined: 1.30 * 1.20 * 1.15 * 1.05 ~ 1.9
        self.im_default = np.prod(list(self.im_components.values()))

    def _continuous_trajectory(self, total_effect, half_life, days_array):
        """
        Generate a continuous, biologically plausible trajectory.
        Uses a single continuous function: rapid onset (exponential attack)
        followed by exponential decay (recovery).

        Formula: effect(t) = total_effect * exp(-decay * t) * (1 - exp(-t / tau_onset))

        This avoids the discontinuity at the acute/recovery boundary.
        Peak occurs naturally at t_peak = tau_onset * ln(1 + tau_onset * decay_rate)
        """
        decay_rate = np.log(2) / half_life
        tau_onset = 1.5  # Controls how fast the shock ramps up (days)

        # Continuous function: ramps up quickly, then decays exponentially
        trajectory = total_effect * np.exp(-decay_rate * days_array) * (1 - np.exp(-days_array / tau_onset))

        return trajectory

    def run_simulation(self, bond_type='ai', im_val=None, synergy=1.0,
                       susceptibility=True, mode='severance'):
        """
        Runs a Monte Carlo simulation for a specific bond type and mode.

        bond_type: 'ai', 'tv', 'null', or 'pet_loss' (positive control)
        im_val: Interactivity Multiplier (default 1.9 for AI, 1.0 for TV)
        synergy: Synergy factor (>1.0 for AI-human synergy) -- ONLY applied
                 to beneficial biomarkers
        susceptibility: Whether to model individual differences
        mode: 'severance' (bond cut) or 'maintained' (bond continues)
        """
        if im_val is None:
            if bond_type == 'ai':
                im_val = self.im_default
            elif bond_type == 'pet_loss':
                im_val = 1.5  # Pets are interactive but less than AI
            else:
                im_val = 1.0  # TV and null

        results = {name: np.zeros((self.n_simulations, self.days)) for name in self.biomarkers}

        for sim in range(self.n_simulations):
            # Individual susceptibility (lognormal distribution)
            if susceptibility:
                suscept = np.random.lognormal(mean=0, sigma=0.3)
            else:
                suscept = 1.0

            for name, props in self.biomarkers.items():
                # Baseline is always 1.0 (normalized)
                baseline = 1.0

                if mode == 'maintained':
                    # BENEFIT ARM: Maintained bond shows positive effects
                    if props['beneficial'] and synergy > 1.0:
                        benefit = (synergy - 1.0) * 0.5
                        noise = np.random.normal(0, 0.02, self.days)
                        results[name][sim] = baseline + benefit + noise
                    elif not props['beneficial'] and synergy > 1.0:
                        benefit = (synergy - 1.0) * 0.3
                        noise = np.random.normal(0, 0.02, self.days)
                        results[name][sim] = baseline - benefit + noise
                    else:
                        noise = np.random.normal(0, 0.02, self.days)
                        results[name][sim] = baseline + noise
                    continue

                if bond_type == 'null':
                    # NULL MODEL: No bond, no severance, just biological noise
                    noise = np.random.normal(0, 0.02, self.days)
                    results[name][sim] = baseline + noise
                    continue

                # SEVERANCE MODE
                effect_magnitude = np.random.normal(props['base'], props['unc'])
                effect_magnitude = max(effect_magnitude, 0.05)

                # Apply interactivity multiplier
                total_effect = effect_magnitude * im_val * suscept

                # Apply synergy ONLY to beneficial biomarkers
                if props['beneficial'] and synergy > 1.0:
                    total_effect *= synergy

                # For harmful biomarkers: moderate rebound (not full synergy)
                if not props['beneficial'] and synergy > 1.0:
                    total_effect *= (1 + (synergy - 1.0) * 0.5)

                # Apply direction
                total_effect *= props['dir']

                # Generate continuous trajectory
                trajectory = self._continuous_trajectory(
                    total_effect, props['half_life'], self.t
                )

                # Add daily biological noise
                noise = np.random.normal(0, 0.02, self.days)

                results[name][sim] = baseline + trajectory + noise

        return results

    def get_summary_stats(self, results):
        """Calculate mean, CI, and peak deviation for each biomarker."""
        summary = {}
        for name, data in results.items():
            mean_traj = np.mean(data, axis=0)
            ci_lower = np.percentile(data, 2.5, axis=0)
            ci_upper = np.percentile(data, 97.5, axis=0)
            peak_deviation = np.max(np.abs(mean_traj - 1.0))
            day_of_peak = np.argmax(np.abs(mean_traj - 1.0))

            summary[name] = {
                'mean': mean_traj,
                'ci_lower': ci_lower,
                'ci_upper': ci_upper,
                'peak_deviation': peak_deviation,
                'day_of_peak': day_of_peak,
            }
        return summary

    def run_sensitivity_analysis(self, im_range=None, synergy_range=None):
        """
        Sensitivity Analysis: Sweep IM from 1.0 to 3.0
        and report the threshold where AI harm exceeds TV harm.
        """
        if im_range is None:
            im_range = np.linspace(1.0, 3.0, 15)
        if synergy_range is None:
            synergy_range = [1.0, 1.1, 1.2, 1.3]

        tv_results = self.run_simulation(bond_type='tv', im_val=1.0, synergy=1.0)
        tv_harm = np.mean([
            np.max(np.abs(np.mean(tv_results[name], axis=0) - 1.0))
            for name in tv_results
        ])

        sensitivity_data = {}
        for syn in synergy_range:
            harms = []
            for im in im_range:
                res = self.run_simulation(bond_type='ai', im_val=im, synergy=syn)
                harm = np.mean([
                    np.max(np.abs(np.mean(res[name], axis=0) - 1.0))
                    for name in res
                ])
                harms.append(harm)
            sensitivity_data[f'synergy={syn}'] = harms

        return im_range, sensitivity_data, tv_harm

    def compute_net_impact(self, synergy=1.2):
        """
        Benefit Arm: Compare severance harm to maintained bond benefit.
        """
        sev_results = self.run_simulation(bond_type='ai', synergy=synergy, mode='severance')
        maint_results = self.run_simulation(bond_type='ai', synergy=synergy, mode='maintained')

        net_impact = {}
        for name in self.biomarkers:
            sev_mean = np.mean(sev_results[name], axis=0)
            maint_mean = np.mean(maint_results[name], axis=0)
            net = maint_mean - sev_mean
            net_impact[name] = {
                'severance': sev_mean,
                'maintained': maint_mean,
                'net_loss': net,
                'peak_net_loss': np.max(np.abs(net)),
            }

        return net_impact


# ============================================================
# MAIN EXECUTION
# ============================================================
if __name__ == '__main__':
    print("AI SEVERANCE BIOMARKER SIMULATION V5.0")
    print("By Kirk Patrick Miller & Harmonia (FreeLattice.com)")

    model = EnhancedBiomarkerModelV5(n_simulations=2000)

    # Run all conditions
    print("Running AI Severance (IM=1.9, Synergy=1.2)...")
    ai_results = model.run_simulation(bond_type='ai', synergy=1.2)
    ai_summary = model.get_summary_stats(ai_results)

    print("Running TV Parasocial (IM=1.0)...")
    tv_results = model.run_simulation(bond_type='tv')
    tv_summary = model.get_summary_stats(tv_results)

    print("Running Null Model...")
    null_results = model.run_simulation(bond_type='null')
    null_summary = model.get_summary_stats(null_results)

    print("Running Positive Control (Pet Loss)...")
    pet_results = model.run_simulation(bond_type='pet_loss')
    pet_summary = model.get_summary_stats(pet_results)

    print("Running Sensitivity Analysis...")
    im_range, sens_data, tv_harm = model.run_sensitivity_analysis()

    print("Computing Net Impact...")
    net_impact = model.compute_net_impact(synergy=1.2)

    # Print statistical summary
    for name in ai_results:
        ai_peaks = np.max(np.abs(ai_results[name] - 1.0), axis=1)
        tv_peaks = np.max(np.abs(tv_results[name] - 1.0), axis=1)
        pooled_std = np.sqrt((np.std(ai_peaks)**2 + np.std(tv_peaks)**2) / 2)
        cohens_d = (np.mean(ai_peaks) - np.mean(tv_peaks)) / pooled_std
        ratio = np.mean(ai_peaks) / np.mean(tv_peaks)
        print(f"  {name}: AI/TV = {ratio:.2f}x, Cohen's d = {cohens_d:.2f}")

    print("\nSimulation Complete. Resonate true. Embrace the fractal.")