Module 2: Neuroscience for EEG Workbook

# Define brain regions and their corresponding electrodes brain_regions = { 'frontal': ['Fp1', 'Fp2', ' ', ' ', 'F3', 'F4', 'F7', 'F8'], 'central': ['C3', 'C4', ' ', ' '], 'parietal': ['P3', 'P4', ' ', ' '], 'temporal': ['T3', 'T4', 'T5', ' '], 'occipital': ['O1', ' '] } # Function to identify brain region from electrode name def get_brain_region(electrode_name): for region, electrodes in brain_regions.items(): if electrode_name in electrodes: return return 'Unknown' # Hemisphere identification def get_hemisphere(electrode_name): if electrode_name[-1] in ['1', '3', '5', '7']: return ' ' elif electrode_name[-1] in ['2', '4', '6', '8']: return ' ' elif electrode_name[-1] == 'z': return ' '

import numpy as np def simulate_neural_dipoles(n_neurons, synchrony_level): """ Simulate how synchronized neural activity creates EEG signals Parameters: - n_neurons: number of neurons - synchrony_level: 0 (random) to 1 (perfectly synchronized) """ # Base frequency of neural oscillation (e.g., 10 Hz for alpha) base_freq = 10 duration = 1 # second fs = 1000 # sampling rate time = np.linspace(0, duration, fs) # Initialize summed activity eeg_signal = np.zeros(len(time)) for i in range(n_neurons): # Random phase offset (reduced by synchrony level) phase_offset = np.random.uniform(0, 2*np.pi) * (1 - ) # Random amplitude (represents dipole strength) amplitude = np.random.normal(1, 0.2) # Individual neuron's contribution neuron_signal = amplitude * np.sin(2 * np.pi * base_freq * time + ) # Add to total EEG eeg_signal += # Scale by number of neurons (spatial summation) eeg_signal = eeg_signal / np.sqrt( ) return time, eeg_signal # What happens to signal amplitude with more synchronized neurons? # Answer:

def detect_brain_rhythms(eeg_data, fs, channel_name): """ Detect dominant brain rhythm in EEG data """ from scipy import signal # Define rhythm bands rhythms = { 'delta': (0.5, 4), 'theta': (4, 8), 'alpha': (8, 13), 'beta': (13, 30), 'gamma': (30, 50) } # Special rhythms by location location_specific = { 'mu': {'freq': (8, 13), 'channels': ['C3', 'C4', ' ']}, 'occipital_alpha': {'freq': (8, 13), 'channels': ['O1', 'O2', ' ']}, 'frontal_theta': {'freq': (4, 8), 'channels': ['Fz', ' ']} } # Compute power spectral density freqs, psd = signal.welch(eeg_data, fs=fs, nperseg=fs*2) # Calculate power in each band band_powers = {} total_power = np.trapz(psd, freqs) for rhythm_name, (low_freq, high_freq) in rhythms.items(): # Find frequency indices idx = np.where((freqs >= low_freq) & (freqs <= ))[0] # Calculate absolute and relative power abs_power = np. (psd[idx], freqs[idx]) rel_power = (abs_power / ) * 100 band_powers[rhythm_name] = { 'absolute': abs_power, 'relative': rel_power } # Find dominant rhythm dominant = max(band_powers.items(), key=lambda x: x[1][' '])[0] # Check for location-specific rhythms if channel_name in location_specific['mu']['channels'] and dominant == 'alpha': # Could be mu rhythm - check for reactivity return 'mu (needs motor task to confirm)' return dominant, band_powers # Alpha reactivity test def test_alpha_reactivity(eeg_eyes_closed, eeg_eyes_open, fs): """Test alpha blocking with eye opening""" # Calculate alpha power for both conditions _, psd_closed = signal.welch(eeg_eyes_closed, fs=fs) _, psd_open = signal.welch(eeg_eyes_open, fs=fs) # Alpha band indices (8-13 Hz) freqs = np.fft.fftfreq(len(psd_closed), 1/fs) alpha_idx = np.where((freqs >= 8) & (freqs <= 13))[0] # Calculate alpha suppression alpha_closed = np.mean(psd_closed[alpha_idx]) alpha_open = np.mean(psd_open[alpha_idx]) suppression_ratio = ( - ) / alpha_closed * 100 return suppression_ratio # Positive value indicates suppression

def calculate_electrode_positions(nasion_to_inion, ear_to_ear): """ Calculate electrode positions based on head measurements Parameters: - nasion_to_inion: distance from nasion to inion (cm) - ear_to_ear: distance between preauricular points (cm) """ # Standard 10-20 percentages positions = {} # Sagittal line (front to back) positions['Fpz'] = nasion_to_inion * # 10% from nasion positions['Fz'] = nasion_to_inion * # 30% from nasion positions['Cz'] = nasion_to_inion * 0.50 # 50% (vertex) positions['Pz'] = nasion_to_inion * # 70% from nasion positions['Oz'] = nasion_to_inion * 0.90 # 90% from nasion # Coronal line (ear to ear) positions['T3'] = ear_to_ear * 0.10 # 10% from left ear positions['C3'] = ear_to_ear * # 30% from left ear positions['Cz_check'] = ear_to_ear * 0.50 # Should match Cz positions['C4'] = ear_to_ear * # 70% from left ear positions['T4'] = ear_to_ear * 0.90 # 90% from left ear return positions # Extended 10-10 system def get_1010_positions(): """Get additional electrodes in 10-10 system""" additional_electrodes = { 'frontal': ['AF3', 'AF4', 'F1', 'F2', 'F5', 'F6'], 'central': ['FC1', 'FC2', 'FC5', 'FC6', 'CP1', 'CP2', 'CP5', 'CP6'], 'parietal': ['P1', 'P2', 'P5', 'P6'], 'temporal': ['FT9', 'FT10', 'TP9', 'TP10'], 'occipital': ['PO3', 'PO4', 'PO7', 'PO8'] } return additional_electrodes # Reference electrode selection reference_types = { 'monopolar': { 'linked_ears': ['A1+A2', 'Good for '], 'Cz': ['Cz', 'Common for studies'], 'average': ['AVG', 'Requires electrodes minimum'] }, 'bipolar': { 'longitudinal': ['Fp1-F3', 'F3-C3', 'C3-P3', 'P3-O1'], 'transverse': ['F7-F3', 'F3-Fz', 'Fz-F4', 'F4-F8'] } }

def create_montage(electrode_data, montage_type): """ Create different montages from monopolar recordings """ if montage_type == 'bipolar_longitudinal': # Define channel pairs pairs = [ ('Fp1', 'F3'), ('F3', 'C3'), ('C3', 'P3'), ('P3', 'O1'), ('Fp2', 'F4'), ('F4', 'C4'), ('C4', 'P4'), ('P4', 'O2'), ('Fp1', 'F7'), ('F7', 'T3'), ('T3', 'T5'), ('T5', 'O1'), ('Fp2', 'F8'), ('F8', 'T4'), ('T4', 'T6'), ('T6', 'O2') ] bipolar_data = {} for ch1, ch2 in pairs: # Bipolar = first electrode minus second bipolar_data[f'{ch1}-{ch2}'] = electrode_data[ch1] - electrode_data[ ] elif montage_type == 'average_reference': # Calculate average of all electrodes all_channels = list(electrode_data.keys()) avg_signal = np.mean([electrode_data[ch] for ch in all_channels], axis= ) avg_ref_data = {} for ch in all_channels: # Subtract average from each channel avg_ref_data[ch] = electrode_data[ch] - elif montage_type == 'laplacian': # Current source density (spatial filter) laplacian_data = {} # Define nearest neighbors for each electrode neighbors = { 'Cz': ['Fz', 'C3', 'C4', 'Pz'], 'C3': ['F3', 'Cz', 'T3', 'P3'], # ... more electrodes } for electrode, neighbor_list in neighbors.items(): if electrode in electrode_data: # Laplacian = electrode - mean of neighbors neighbor_mean = np.mean([electrode_data[n] for n in neighbor_list if n in electrode_data], axis=0) laplacian_data[electrode] = electrode_data[electrode] - return bipolar_data if montage_type.startswith('bipolar') else avg_ref_data

def detect_sleep_features(eeg_data, fs, stage='N2'): """ Detect characteristic features of different sleep stages """ if stage == 'N2': # Detect sleep spindles spindles = detect_spindles(eeg_data, fs) # Detect K-complexes k_complexes = detect_k_complexes(eeg_data, fs) return {'spindles': spindles, 'k_complexes': k_complexes} elif stage == 'N3': # Detect slow waves slow_waves = detect_slow_waves(eeg_data, fs) return {'slow_waves': slow_waves} def detect_spindles(eeg_data, fs): """ Detect sleep spindles (12-14 Hz, 0.5-2 sec duration) """ from scipy import signal # Bandpass filter for spindle frequency spindle_band = ( , ) # Hz nyquist = fs / 2 b, a = signal.butter(4, [spindle_band[0]/nyquist, spindle_band[1]/nyquist], 'band') # Filter the signal spindle_filtered = signal.filtfilt(b, a, eeg_data) # Calculate envelope using Hilbert transform analytic_signal = signal.hilbert(spindle_filtered) envelope = np.abs( ) # Threshold for spindle detection (e.g., 2 SD above mean) threshold = np.mean(envelope) + * np.std(envelope) # Find segments above threshold above_threshold = envelope > threshold # Identify spindle events (duration criteria) min_duration = * fs # 0.5 seconds max_duration = * fs # 2 seconds spindle_events = [] in_spindle = False start_idx = 0 for i, val in enumerate(above_threshold): if val and not in_spindle: in_spindle = True start_idx = i elif not val and in_spindle: in_spindle = False duration = i - start_idx if min_duration <= duration <= max_duration: spindle_events.append((start_idx/fs, i/fs)) # Convert to seconds return spindle_events def detect_k_complexes(eeg_data, fs): """ Detect K-complexes (sharp negative then positive wave) """ # Low-pass filter to enhance slow components cutoff = 4 # Hz b, a = signal.butter(4, cutoff/(fs/2), 'low') filtered = signal.filtfilt(b, a, eeg_data) # Find large amplitude deflections # K-complex: negative peak followed by positive within 0.5-1 sec # Calculate derivative to find sharp changes derivative = np.diff(filtered) # Find negative peaks neg_peaks, _ = signal.find_peaks(-filtered, height= , # Minimum amplitude (μV) distance=int(1.5*fs)) # Minimum 1.5 sec apart k_complexes = [] for neg_peak in neg_peaks: # Look for positive deflection within 0.5-1 second window_start = neg_peak window_end = min(neg_peak + int(1*fs), len(filtered)) window_data = filtered[window_start:window_end] if len(window_data) > 0: pos_peak = np.argmax(window_data) # Check if positive peak is significant if window_data[pos_peak] > : # μV threshold k_complexes.append(neg_peak/fs) # Convert to seconds return k_complexes

class EEGRhythmAnalyzer: def __init__(self, fs=256): self.fs = fs self.rhythms = { 'delta': {'range': (0.5, 4), 'amplitude': (20, 200)}, 'theta': {'range': (4, 8), 'amplitude': (10, 100)}, 'alpha': {'range': (8, 13), 'amplitude': (20, 60)}, 'beta': {'range': (13, 30), 'amplitude': (5, 30)}, 'gamma': {'range': (30, 100), 'amplitude': (1, 5)} } def analyze_channel(self, eeg_data, channel_name): """ Comprehensive rhythm analysis for a single channel """ results = { 'channel': channel_name, 'dominant_rhythm': None, 'rhythm_powers': {}, 'special_features': [] } # Calculate PSD freqs, psd = signal.welch(eeg_data, fs=self.fs, nperseg=self.fs*2) # Analyze each rhythm max_power = 0 dominant = None for rhythm_name, params in self.rhythms.items(): freq_range = params['range'] # Extract band band_idx = np.where((freqs >= freq_range[0]) & (freqs <= freq_range[1]))[0] # Calculate metrics band_power = np.trapz(psd[band_idx], freqs[band_idx]) peak_freq = freqs[band_idx[np.argmax(psd[band_idx])]] peak_amplitude = np.sqrt(np.max(psd[band_idx]) * 2) # Convert to amplitude results['rhythm_powers'][rhythm_name] = { 'power': band_power, 'peak_frequency': peak_freq, 'peak_amplitude': peak_amplitude, 'in_normal_range': params['amplitude'][0] <= peak_amplitude <= params['amplitude'][1] } if band_power > max_power: max_power = band_power dominant = rhythm_name results['dominant_rhythm'] = dominant # Check for special features based on location if channel_name in ['O1', 'O2'] and dominant == 'alpha': results['special_features'].append('Posterior dominant rhythm') elif channel_name in ['C3', 'C4'] and dominant == 'alpha': # Could be mu rhythm - check characteristics if self._check_mu_characteristics(eeg_data): results['special_features'].append('Possible mu rhythm') elif channel_name == 'Fz' and results['rhythm_powers']['theta']['power'] > \ results['rhythm_powers']['alpha']['power']: results['special_features'].append('Frontal midline theta') return results def _check_mu_characteristics(self, eeg_data): """ Check if alpha-range activity has mu rhythm characteristics """ # Mu rhythm has arch-shaped morphology # Check for non-sinusoidal shape using spectral analysis freqs, psd = signal.periodogram(eeg_data, fs=self.fs) # Find 10 Hz peak peak_10hz_idx = np.argmin(np.abs(freqs - 10)) # Check for harmonics (mu has strong 20 Hz harmonic) peak_20hz_idx = np.argmin(np.abs(freqs - 20)) harmonic_ratio = psd[peak_20hz_idx] / psd[peak_10hz_idx] # Mu rhythm typically has stronger harmonics than alpha return harmonic_ratio > # Threshold value def compare_conditions(self, eeg_rest, eeg_task, channel): """ Compare rhythms between two conditions (e.g., rest vs task) """ rest_analysis = self.analyze_channel(eeg_rest, channel) task_analysis = self.analyze_channel(eeg_task, channel) comparison = { 'channel': channel, 'changes': {} } for rhythm in self.rhythms.keys(): rest_power = rest_analysis['rhythm_powers'][rhythm]['power'] task_power = task_analysis['rhythm_powers'][rhythm]['power'] change_percent = ((task_power - rest_power) / rest_power) * 100 comparison['changes'][rhythm] = { 'rest_power': rest_power, 'task_power': task_power, 'change_percent': change_percent, 'significant': abs(change_percent) > # 20% threshold } # Identify reactivity patterns if channel in ['O1', 'O2'] and comparison['changes']['alpha']['change_percent'] < -30: comparison['reactivity'] = 'Alpha blocking detected' elif channel in ['C3', 'C4'] and comparison['changes']['alpha']['change_percent'] < -30: comparison['reactivity'] = 'Mu suppression detected' return comparison