synthetic-sciences Pharmacology Wet-Lab: Experimental Data Analysis
Overview
Pharmacology Wet-Lab provides computational tools for analyzing data from pharmacology experiments. This skill covers western blot densitometry and quantification, xenograft tumor growth inhibition analysis, pharmaceutical stability modeling using Arrhenius kinetics, radiolabeled antibody biodistribution calculations, MIRD-based dosimetry, adverse event grading against CTCAE criteria, and dose-response curve fitting for IC50/EC50 determination.
When to Use This Skill
- Quantifying protein expression from western blot images
- Analyzing xenograft tumor growth data and calculating TGI%
- Predicting pharmaceutical shelf life from accelerated stability data
- Processing radiolabeled antibody biodistribution data (%ID/g)
- Estimating absorbed radiation doses (MIRD dosimetry)
- Grading adverse events against CTCAE or VCOG-CTCAE scales
- Fitting dose-response curves for IC50/EC50 determination
- Calculating combination indices (Chou-Talalay method)
Related Skills: For drug database queries use chembl-database or fda-database. For molecular docking use diffdock. For survival analysis use scikit-survival.
Installation
uv pip install opencv-python scipy pandas numpy matplotlib lifelines
Quick Start
import numpy as np
from scipy.optimize import curve_fit
# 4-Parameter Logistic for dose-response (IC50)
def four_pl(x, bottom, top, ic50, hill):
return bottom + (top - bottom) / (1 + (x / ic50) ** hill)
concentrations = np.array([0.001, 0.01, 0.1, 1, 10, 100]) # uM
viability = np.array([98, 95, 82, 45, 12, 3]) # % viability
popt, pcov = curve_fit(four_pl, concentrations, viability,
p0=[0, 100, 1, 1], maxfev=10000)
print(f"IC50: {popt[2]:.3f} uM")
print(f"Hill coefficient: {popt[3]:.2f}")
Core Capabilities
1. Western Blot Densitometry
Quantify protein bands from western blot images.
import cv2
import numpy as np
def quantify_western_blot(image_path, n_lanes, band_height=50):
"""Quantify western blot band intensities.
Args:
image_path: path to blot image
n_lanes: number of lanes
band_height: expected band height in pixels
"""
image = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
if image is None:
raise FileNotFoundError(f"Cannot load {image_path}")
# Invert (bands are dark on light background)
inverted = 255 - image
h, w = inverted.shape
# Divide into lanes
lane_width = w // n_lanes
lane_intensities = []
for i in range(n_lanes):
x_start = i * lane_width + int(lane_width * 0.1)
x_end = (i + 1) * lane_width - int(lane_width * 0.1)
lane = inverted[:, x_start:x_end]
# Find band (peak in vertical intensity profile)
profile = lane.mean(axis=1)
# Background subtraction (rolling ball approximation)
from scipy.ndimage import minimum_filter1d
background = minimum_filter1d(profile, size=100)
corrected = profile - background
corrected = np.clip(corrected, 0, None)
# Band detection
from scipy.signal import find_peaks
peaks, props = find_peaks(corrected, height=corrected.max()*0.1,
distance=band_height)
# Integrate band intensity (area under curve)
total_intensity = 0
for peak in peaks:
start = max(0, peak - band_height // 2)
end = min(len(corrected), peak + band_height // 2)
band_area = corrected[start:end].sum()
total_intensity += band_area
lane_intensities.append({
'lane': i + 1,
'raw_intensity': total_intensity,
'n_bands': len(peaks),
'peak_positions': list(peaks)
})
# Normalize to loading control (first lane or specified)
import pandas as pd
df = pd.DataFrame(lane_intensities)
control_intensity = df.iloc[0]['raw_intensity']
df['normalized'] = df['raw_intensity'] / control_intensity
df['fold_change'] = df['normalized']
print("Lane intensities:")
print(df[['lane', 'raw_intensity', 'normalized', 'fold_change']])
return df
def calculate_fold_change(target_intensities, loading_control_intensities):
"""Calculate normalized fold change with loading control.
Args:
target_intensities: list of target protein band intensities
loading_control_intensities: list of loading control (e.g., actin) intensities
"""
target = np.array(target_intensities)
control = np.array(loading_control_intensities)
# Normalize target to loading control
normalized = target / control
# Fold change relative to first sample
fold_change = normalized / normalized[0]
return fold_change
2. Xenograft Tumor Growth Analysis
Analyze in vivo tumor growth and treatment efficacy.
import numpy as np
import pandas as pd
from scipy import stats
def calculate_tgi(tumor_volumes, groups, time_points, control_group='Vehicle'):
"""Calculate Tumor Growth Inhibition (TGI%).
TGI% = (1 - (Tt - T0) / (Ct - C0)) * 100
where T=treatment, C=control, t=final, 0=initial
Args:
tumor_volumes: DataFrame with columns [animal_id, group, time, volume]
control_group: name of control group
"""
results = []
final_time = max(time_points)
# Control group growth
ctrl = tumor_volumes[tumor_volumes['group'] == control_group]
ctrl_initial = ctrl[ctrl['time'] == 0]['volume'].mean()
ctrl_final = ctrl[ctrl['time'] == final_time]['volume'].mean()
ctrl_growth = ctrl_final - ctrl_initial
for group in tumor_volumes['group'].unique():
if group == control_group:
continue
treat = tumor_volumes[tumor_volumes['group'] == group]
treat_initial = treat[treat['time'] == 0]['volume'].mean()
treat_final = treat[treat['time'] == final_time]['volume'].mean()
treat_growth = treat_final - treat_initial
tgi = (1 - treat_growth / ctrl_growth) * 100
# Statistical comparison at final timepoint
ctrl_final_vals = ctrl[ctrl['time'] == final_time]['volume'].values
treat_final_vals = treat[treat['time'] == final_time]['volume'].values
stat, pval = stats.mannwhitneyu(ctrl_final_vals, treat_final_vals,
alternative='two-sided')
results.append({
'group': group,
'TGI%': tgi,
'mean_volume_final': treat_final,
'control_volume_final': ctrl_final,
'p_value': pval,
'n_animals': len(treat['animal_id'].unique())
})
df = pd.DataFrame(results)
print("Tumor Growth Inhibition:")
print(df[['group', 'TGI%', 'mean_volume_final', 'p_value']])
return df
def growth_delay(tumor_volumes, groups, doubling_volume):
"""Calculate tumor growth delay (time to reach doubling volume).
Args:
doubling_volume: target volume for comparison (e.g., 2x initial)
"""
results = []
for group in tumor_volumes['group'].unique():
grp = tumor_volumes[tumor_volumes['group'] == group]
times_to_double = []
for animal in grp['animal_id'].unique():
animal_data = grp[grp['animal_id'] == animal].sort_values('time')
reached = animal_data[animal_data['volume'] >= doubling_volume]
if len(reached) > 0:
times_to_double.append(reached.iloc[0]['time'])
if times_to_double:
results.append({
'group': group,
'median_time_to_double': np.median(times_to_double),
'mean_time_to_double': np.mean(times_to_double),
'n_reached': len(times_to_double)
})
return pd.DataFrame(results)
3. Pharmaceutical Stability
Predict shelf life from accelerated stability data.
import numpy as np
from scipy.optimize import curve_fit
def arrhenius_stability(temperatures_C, rate_constants):
"""Fit Arrhenius equation to predict shelf life.
k = A * exp(-Ea / RT)
Args:
temperatures_C: storage temperatures in Celsius
rate_constants: degradation rate constants at each temperature
"""
T_kelvin = np.array(temperatures_C) + 273.15
k = np.array(rate_constants)
# Linearize: ln(k) = ln(A) - Ea/(R*T)
ln_k = np.log(k)
inv_T = 1 / T_kelvin
# Linear regression
slope, intercept, r_value, p_value, std_err = stats.linregress(inv_T, ln_k)
R = 8.314 # J/(mol*K)
Ea = -slope * R / 1000 # kJ/mol
A = np.exp(intercept)
print(f"Activation energy: {Ea:.1f} kJ/mol")
print(f"Pre-exponential factor: {A:.4e}")
print(f"R²: {r_value**2:.4f}")
# Predict rate at storage temperature (25°C)
T_storage = 25 + 273.15
k_25 = A * np.exp(-Ea * 1000 / (R * T_storage))
# Shelf life (time to 10% degradation, first-order kinetics)
# C/C0 = exp(-k*t) → t = -ln(0.9) / k
shelf_life_days = -np.log(0.9) / k_25
shelf_life_months = shelf_life_days / 30.44
print(f"\nPredicted at 25°C:")
print(f" Rate constant: {k_25:.6f} day⁻¹")
print(f" Shelf life (to 90%): {shelf_life_months:.0f} months ({shelf_life_days:.0f} days)")
return {
'Ea_kJ_mol': Ea,
'A': A,
'k_25C': k_25,
'shelf_life_months': shelf_life_months,
'shelf_life_days': shelf_life_days
}
def fit_degradation_kinetics(time_points, concentrations, order=1):
"""Fit degradation kinetics (zero, first, or second order)."""
time = np.array(time_points)
conc = np.array(concentrations)
if order == 0:
# C = C0 - k*t
slope, intercept, r, p, se = stats.linregress(time, conc)
k = -slope
r_sq = r**2
elif order == 1:
# ln(C) = ln(C0) - k*t
slope, intercept, r, p, se = stats.linregress(time, np.log(conc))
k = -slope
r_sq = r**2
elif order == 2:
# 1/C = 1/C0 + k*t
slope, intercept, r, p, se = stats.linregress(time, 1/conc)
k = slope
r_sq = r**2
print(f"Order {order}: k = {k:.6f}, R² = {r_sq:.4f}")
return {'k': k, 'r_squared': r_sq, 'order': order}
4. Radiolabeled Antibody Biodistribution
Process biodistribution data from radiolabeled compounds.
import numpy as np
import pandas as pd
def calculate_biodistribution(tissue_counts, tissue_weights_g, injection_dose,
decay_correction_factor=1.0):
"""Calculate %ID/g from tissue radioactivity counts.
Args:
tissue_counts: dict of tissue_name -> counts per minute (CPM)
tissue_weights_g: dict of tissue_name -> weight in grams
injection_dose: total injected dose in CPM
decay_correction_factor: factor to correct for radioactive decay
"""
results = []
for tissue in tissue_counts:
cpm = tissue_counts[tissue] * decay_correction_factor
weight = tissue_weights_g[tissue]
id_per_g = (cpm / weight) / injection_dose * 100 # %ID/g
results.append({
'tissue': tissue,
'CPM': tissue_counts[tissue],
'weight_g': weight,
'percent_ID_per_g': id_per_g
})
df = pd.DataFrame(results).sort_values('percent_ID_per_g', ascending=False)
# Tumor-to-blood ratio (if both present)
if 'tumor' in tissue_counts and 'blood' in tissue_counts:
tumor_uptake = df[df['tissue'] == 'tumor']['percent_ID_per_g'].values[0]
blood_uptake = df[df['tissue'] == 'blood']['percent_ID_per_g'].values[0]
ratio = tumor_uptake / blood_uptake if blood_uptake > 0 else float('inf')
print(f"Tumor-to-blood ratio: {ratio:.2f}")
print(df[['tissue', 'percent_ID_per_g']].to_string(index=False))
return df
5. MIRD Dosimetry
Estimate absorbed radiation dose using MIRD formalism.
import numpy as np
from scipy.integrate import trapezoid
def mird_dosimetry(time_points_h, activity_values_MBq, s_values, organ_masses_kg):
"""Estimate absorbed dose using MIRD formalism.
Dose = cumulative_activity * S_value
Args:
time_points_h: time points in hours
activity_values_MBq: dict of organ -> list of activity values at each time
s_values: dict of (source, target) -> S-value in mGy/(MBq*h)
organ_masses_kg: dict of organ -> mass in kg
"""
results = []
time = np.array(time_points_h)
for organ, activities in activity_values_MBq.items():
act = np.array(activities)
# Time-integrated activity (cumulative activity, A_tilde)
cumulative_activity = trapezoid(act, time) # MBq*h
# Self-dose
s_self = s_values.get((organ, organ), 0)
dose_self = cumulative_activity * s_self # mGy
# Cross-dose from other organs
dose_cross = 0
for source_organ, source_activities in activity_values_MBq.items():
if source_organ != organ:
s_cross = s_values.get((source_organ, organ), 0)
source_cumulative = trapezoid(np.array(source_activities), time)
dose_cross += source_cumulative * s_cross
total_dose = dose_self + dose_cross
results.append({
'organ': organ,
'cumulative_activity_MBq_h': cumulative_activity,
'self_dose_mGy': dose_self,
'cross_dose_mGy': dose_cross,
'total_dose_mGy': total_dose
})
import pandas as pd
df = pd.DataFrame(results)
print(df[['organ', 'cumulative_activity_MBq_h', 'total_dose_mGy']])
return df
6. Adverse Event Grading
Grade adverse events against CTCAE v5.0 criteria.
import pandas as pd
# CTCAE v5.0 grading examples (subset)
CTCAE_CRITERIA = {
'neutrophil_count': {
1: lambda x: 1500 <= x < 2000, # Grade 1: LLN - 1500/mm3
2: lambda x: 1000 <= x < 1500, # Grade 2: 1000-1500
3: lambda x: 500 <= x < 1000, # Grade 3: 500-1000
4: lambda x: x < 500, # Grade 4: <500
},
'platelet_count': {
1: lambda x: 75000 <= x < 150000,
2: lambda x: 50000 <= x < 75000,
3: lambda x: 25000 <= x < 50000,
4: lambda x: x < 25000,
},
'hemoglobin': {
1: lambda x: 10 <= x < 12, # g/dL
2: lambda x: 8 <= x < 10,
3: lambda x: x < 8,
4: lambda x: False, # Life-threatening
},
'ALT': {
1: lambda x: 1 <= x/40 < 3, # x ULN (ULN=40)
2: lambda x: 3 <= x/40 < 5,
3: lambda x: 5 <= x/40 < 20,
4: lambda x: x/40 >= 20,
}
}
def grade_adverse_event(parameter, value):
"""Grade an adverse event by CTCAE v5.0."""
if parameter not in CTCAE_CRITERIA:
return {'grade': None, 'error': f'Unknown parameter: {parameter}'}
criteria = CTCAE_CRITERIA[parameter]
for grade in [4, 3, 2, 1]: # Check highest grade first
if criteria[grade](value):
return {'parameter': parameter, 'value': value, 'grade': grade}
return {'parameter': parameter, 'value': value, 'grade': 0}
def assess_dlt(adverse_events, dlt_grades={'hematologic': 4, 'non_hematologic': 3}):
"""Assess dose-limiting toxicity from adverse event list."""
dlts = []
for ae in adverse_events:
if ae['grade'] >= dlt_grades.get('hematologic', 4):
dlts.append(ae)
return dlts
7. Dose-Response Curves
Fit dose-response data for IC50/EC50 determination.
import numpy as np
from scipy.optimize import curve_fit
def four_pl(x, bottom, top, ec50, hill):
"""4-Parameter Logistic dose-response."""
return bottom + (top - bottom) / (1 + (x / ec50) ** hill)
def fit_dose_response(concentrations, responses, response_type='inhibition'):
"""Fit dose-response curve and calculate IC50/EC50.
Args:
concentrations: drug concentrations
responses: measured responses (viability %, activity, etc.)
response_type: 'inhibition' (response decreases) or 'activation'
"""
conc = np.array(concentrations)
resp = np.array(responses)
if response_type == 'inhibition':
p0 = [resp.min(), resp.max(), np.median(conc), 1.0]
else:
p0 = [resp.min(), resp.max(), np.median(conc), 1.0]
bounds = ([0, 0, 0, 0.1], [200, 200, conc.max() * 10, 10])
popt, pcov = curve_fit(four_pl, conc, resp, p0=p0, bounds=bounds, maxfev=10000)
perr = np.sqrt(np.diag(pcov))
bottom, top, ec50, hill = popt
# R-squared
predicted = four_pl(conc, *popt)
ss_res = np.sum((resp - predicted) ** 2)
ss_tot = np.sum((resp - np.mean(resp)) ** 2)
r_squared = 1 - ss_res / ss_tot
label = 'IC50' if response_type == 'inhibition' else 'EC50'
print(f"{label}: {ec50:.4f} ± {perr[2]:.4f}")
print(f"Hill coefficient: {hill:.2f}")
print(f"Bottom: {bottom:.1f}%, Top: {top:.1f}%")
print(f"R²: {r_squared:.4f}")
return {label: ec50, 'hill': hill, 'bottom': bottom, 'top': top,
'r_squared': r_squared, 'params': popt}
def combination_index(drug1_ic50, drug2_ic50, combo_ic50_d1, combo_ic50_d2):
"""Calculate Chou-Talalay Combination Index.
CI < 1: synergistic, CI = 1: additive, CI > 1: antagonistic
"""
ci = (combo_ic50_d1 / drug1_ic50) + (combo_ic50_d2 / drug2_ic50)
interpretation = 'synergistic' if ci < 0.9 else ('antagonistic' if ci > 1.1 else 'additive')
print(f"CI = {ci:.3f} ({interpretation})")
return {'CI': ci, 'interpretation': interpretation}
Typical Workflows
Workflow 1: Quantify Western Blot Bands and Calculate Fold Change
# Target protein bands
target = [1500, 3200, 4800, 2100] # Arbitrary units per lane
actin = [2000, 1950, 2100, 1900] # Loading control
fc = calculate_fold_change(target, actin)
for i, f in enumerate(fc):
print(f"Lane {i+1}: fold change = {f:.2f}")
Workflow 2: Analyze Xenograft Tumor Growth and Calculate TGI%
import pandas as pd
data = pd.DataFrame({
'animal_id': ['V1','V1','V2','V2','T1','T1','T2','T2'],
'group': ['Vehicle','Vehicle','Vehicle','Vehicle','Drug','Drug','Drug','Drug'],
'time': [0, 14, 0, 14, 0, 14, 0, 14],
'volume': [100, 800, 120, 750, 110, 350, 105, 380]
})
tgi = calculate_tgi(data, data['group'].unique(), [0, 14])
Workflow 3: Predict Shelf Life from Accelerated Stability Data
from scipy import stats
temps = [40, 50, 60, 70] # Celsius
rates = [0.001, 0.003, 0.01, 0.03] # day^-1
result = arrhenius_stability(temps, rates)
print(f"Predicted shelf life at 25°C: {result['shelf_life_months']:.0f} months")
Best Practices
- Western blots — always normalize to loading control (beta-actin, GAPDH, total protein); report fold change, not raw intensity
- Xenograft TGI — require minimum 8-10 animals per group; report both TGI% and statistical significance; use caliper measurements consistently
- Stability studies — use at least 3 temperatures for Arrhenius fitting; verify linearity of ln(k) vs 1/T; ICH guidelines require 25°C/60%RH long-term
- Dose-response — use at least 6 concentrations spanning 3 log units; include 0 and saturating dose; fit with Hill slope unconstrained
- Combination index — measure IC50 of each drug alone and in fixed-ratio combination; CI only valid at Fa (fraction affected) near 0.5
- Adverse events — use CTCAE v5.0 for human clinical trials; VCOG-CTCAE for veterinary studies
Troubleshooting
Problem: Dose-response fit gives unreasonable IC50 Solution: Ensure concentrations span the IC50 (responses should range from ~10% to ~90% effect). Add more points around the inflection. Check if response truly follows sigmoidal pattern.
Problem: Western blot quantification inconsistent Solution: Ensure uniform exposure across lanes. Use ECL with linear dynamic range. Avoid saturated bands. Use total protein stain (Ponceau S) instead of housekeeping gene for normalization.
Problem: Arrhenius fit has low R² Solution: Verify degradation follows the assumed kinetic order. Check for multiple degradation mechanisms at different temperatures. Exclude temperatures where phase transitions occur.
Problem: TGI% is negative Solution: Treatment group grew faster than control. Verify group assignments. Check for handling errors. This can occur with growth factors or hormones.