# R
# R: Cell Biology Timeline Visualization
# Mapping the major discoveries that shaped our understanding of cells

# Load required libraries
if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
if (!requireNamespace("ggrepel", quietly = TRUE)) install.packages("ggrepel")

library(ggplot2)
library(ggrepel)

# Historical milestones in cell biology
milestones <- data.frame(
  year = c(1665, 1674, 1838, 1839, 1855, 1931, 1953, 1973, 1998, 2020),
  scientist = c("Hooke", "Leeuwenhoek", "Schleiden", "Schwann", "Virchow", 
                "Ruska", "Watson & Crick", "Boyer", "Yamanaka", "CRISPR 2.0"),
  discovery = c("First cells observed", "Living microorganisms", "Plant cell theory", 
                "Animal cell theory", "Cell division", "Electron microscopy", 
                "DNA structure", "Genetic engineering", "iPSCs", "Base editing"),
  impact = c(8, 9, 10, 10, 9, 8, 10, 9, 8, 7),
  category = c("Microscopy", "Microscopy", "Theory", "Theory", "Theory", 
               "Technology", "Molecular", "Biotech", "Stem Cells", "Gene Editing")
)

# Create timeline visualization
timeline_plot <- ggplot(milestones, aes(x = year, y = impact, color = category)) +
  geom_point(size = 4, alpha = 0.8) +
  geom_line(aes(group = 1), alpha = 0.3, linewidth = 1) +  # use linewidth (not size)
  geom_text_repel(aes(label = paste(scientist, "\n", discovery)), 
                  size = 3, box.padding = 0.5, 
                  point.padding = 0.3, segment.alpha = 0.6) +
  scale_x_continuous(breaks = seq(1650, 2020, 50), limits = c(1650, 2030)) +
  scale_y_continuous(breaks = 1:10, labels = c("1","2","3","4","5","6","7","8","9","Revolutionary")) +
  labs(title = "Cell Biology: 355 Years of Discovery",
       subtitle = "From Hooke's Cork Cells to CRISPR Gene Editing",
       x = "Year of Discovery",
       y = "Impact on Cell Biology",
       color = "Research Area") +
  theme_minimal() +
  theme(
    plot.title = element_text(size = 16, face = "bold"),
    plot.subtitle = element_text(size = 12, face = "italic"),
    legend.position = "bottom"
  )

print(timeline_plot)

cat("\n=== Cell Biology Discovery Timeline ===\n")
cat("Time span:", max(milestones$year) - min(milestones$year), "years\n")
cat("Average impact score:", round(mean(milestones$impact), 2), "\n")
cat("Most revolutionary periods:\n")

revolutionary_periods <- milestones[milestones$impact >= 9, ]
for (i in 1:nrow(revolutionary_periods)) {
  cat("  ", revolutionary_periods$year[i], ":", revolutionary_periods$discovery[i], "\n")
}

# R
# R: Cellular Diversity Analysis
# Comparing characteristics across the spectrum of life

# Only load what's needed, quietly
suppressPackageStartupMessages({
  if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
  if (!requireNamespace("ggrepel", quietly = TRUE)) install.packages("ggrepel")
  library(ggplot2)
  library(ggrepel)
})

# Comprehensive cellular diversity dataset
cellular_diversity <- data.frame(
  organism = c("Mycoplasma", "E. coli", "Yeast", "Human RBC", "Neuron",
               "Plant cell", "Amoeba", "Ostrich egg"),
  type = c("Prokaryote", "Prokaryote", "Eukaryote", "Eukaryote", "Eukaryote",
           "Eukaryote", "Eukaryote", "Eukaryote"),
  diameter_um = c(0.2, 2, 5, 8, 100, 30, 500, 170000),
  volume_um3 = c(0.004, 4, 65, 268, 523000, 14000, 65000000, 2.57e12),
  genome_mb = c(0.16, 4.6, 12.1, 0, 3200, 125, 670, 3200),
  complexity = c(1, 2, 4, 3, 9, 6, 5, 8),
  lifespan_days = c(1, 0.04, 7, 120, 36500, 3650, 30, 1)
)

# Simple 2D visualization with ggplot2 viridis scale (no plotly needed)
diversity_plot <- ggplot(cellular_diversity, aes(x = log10(diameter_um), y = log10(volume_um3),
                                                 color = type, size = complexity)) +
  geom_point(alpha = 0.8) +
  geom_text_repel(aes(label = organism), size = 3) +
  scale_color_viridis_d() +
  labs(title = "Cellular Diversity Across Life",
       x = "log₁₀ Diameter (μm)",
       y = "log₁₀ Volume (μm³)",
       color = "Cell Type",
       size = "Complexity") +
  theme_minimal()

print(diversity_plot)

# Statistical analysis of cellular diversity
cat("\n=== Cellular Diversity Statistics ===\n")
cat("Size range (diameter):", min(cellular_diversity$diameter_um), "to",
    max(cellular_diversity$diameter_um), "μm\n")
cat("Size ratio (largest/smallest):",
    round(max(cellular_diversity$diameter_um)/min(cellular_diversity$diameter_um)), "\n")
cat("Volume range:", formatC(min(cellular_diversity$volume_um3), format = "e", digits = 2),
    "to", formatC(max(cellular_diversity$volume_um3), format = "e", digits = 2), "μm³\n")
cat("Genome size range:",
    min(cellular_diversity$genome_mb[cellular_diversity$genome_mb > 0]),
    "to", max(cellular_diversity$genome_mb), "Mb\n")

# Cell type distribution
cell_types <- table(cellular_diversity$type)
cat("\nCell type distribution:\n")
print(cell_types)

# R: Surface Area to Volume Analysis
# Understanding the fundamental constraint on cell size

# Function to calculate SA:V ratio for spherical cells
sa_to_vol_ratio <- function(radius) {
  surface_area <- 4 * pi * radius^2
  volume <- (4/3) * pi * radius^3
  ratio <- surface_area / volume
  return(ratio)
}

# Cell radius range (micrometers)
radii <- seq(0.1, 50, by = 0.5)
ratios <- sapply(radii, sa_to_vol_ratio)

# Create data frame
sa_data <- data.frame(radius = radii, sa_vol_ratio = ratios)

# Add cell type categories
sa_data$cell_type <- ifelse(sa_data$radius < 1, "Bacteria",
                      ifelse(sa_data$radius < 10, "Small Eukaryote",
                      ifelse(sa_data$radius < 25, "Large Eukaryote", "Giant Cell")))

# Plot SA:V ratio vs cell radius
p1 <- ggplot(sa_data, aes(x = radius, y = sa_vol_ratio, color = cell_type)) +
    geom_line(size = 1.5) +
    scale_y_log10() +
    labs(title = "Surface Area to Volume Ratio vs Cell Radius",
         x = "Cell Radius (μm)",
         y = "SA:V Ratio (μm⁻¹, log scale)",
         color = "Cell Type") +
    theme_minimal() +
    geom_hline(yintercept = 1, linetype = "dashed", color = "red", alpha = 0.7) +
    annotate("text", x = 40, y = 1.2, label = "Critical efficiency threshold", color = "red")

print(p1)

# Real cell examples
real_cells <- data.frame(
    cell_type = c("E. coli", "Yeast", "Red blood cell", "Muscle cell", "Ostrich egg"),
    radius_um = c(1.0, 2.5, 4.0, 50, 85000),
    sa_vol_ratio = sapply(c(1.0, 2.5, 4.0, 50, 85000), sa_to_vol_ratio)
)

# Plot real cell examples
p2 <- ggplot(real_cells, aes(x = reorder(cell_type, radius_um), y = sa_vol_ratio)) +
    geom_col(fill = "steelblue", alpha = 0.7) +
    scale_y_log10() +
    labs(title = "SA:V Ratios of Real Cells",
         x = "Cell Type (ordered by size)",
         y = "SA:V Ratio (μm⁻¹, log scale)") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    geom_text(aes(label = round(sa_vol_ratio, 3)), vjust = -0.2, size = 3)

print(p2)

cat("\n=== Surface Area to Volume Analysis ===\n")
cat("Why cells stay small:\n")
cat("- Small cells have higher SA:V ratios\n")
cat("- Better nutrient uptake and waste removal\n")
cat("- More efficient cellular processes\n\n")

cat("Real cell SA:V ratios:\n")
for(i in 1:nrow(real_cells)) {
    cat(sprintf("  %s (r=%.1f μm): SA:V = %.3f\n", 
                real_cells$cell_type[i], 
                real_cells$radius_um[i], 
                real_cells$sa_vol_ratio[i]))
}

# Efficiency categories
high_eff <- sum(real_cells$sa_vol_ratio > 1)
low_eff <- sum(real_cells$sa_vol_ratio <= 1)
cat("\nEfficiency distribution:\n")
cat("  High efficiency (SA:V > 1):", high_eff, "cell types\n")
cat("  Low efficiency (SA:V ≤ 1):", low_eff, "cell types\n")

# R: Prokaryotic Growth Analysis
# Modeling bacterial population dynamics

# Exponential growth model
exponential_growth <- function(N0, r, t) {
    N0 * exp(r * t)
}

# Logistic growth model
logistic_growth <- function(N0, r, K, t) {
    (K * N0 * exp(r * t)) / (K + N0 * (exp(r * t) - 1))
}

# Growth parameters for different bacteria
bacteria_data <- data.frame(
    species = c("E. coli", "B. subtilis", "S. aureus", "M. tuberculosis"),
    doubling_time_min = c(20, 30, 28, 720),  # minutes
    growth_rate = log(2) / c(20, 30, 28, 720) * 60,  # per hour
    optimal_temp = c(37, 37, 37, 37),
    environment = c("Intestinal", "Soil", "Pathogenic", "Pathogenic")
)

# Simulation parameters
time_hours <- seq(0, 24, 0.5)
N0 <- 100  # Initial population
K <- 1e9   # Carrying capacity

# Generate growth curves for first 3 species (M. tuberculosis is too slow)
growth_curves <- data.frame()
for(i in 1:3) {
    species <- bacteria_data$species[i]
    r <- bacteria_data$growth_rate[i]
    
    exp_pop <- exponential_growth(N0, r, time_hours)
    log_pop <- logistic_growth(N0, r, K, time_hours)
    
    temp_data <- data.frame(
        time = rep(time_hours, 2),
        population = c(exp_pop, log_pop),
        model = rep(c("Exponential", "Logistic"), each = length(time_hours)),
        species = species
    )
    
    growth_curves <- rbind(growth_curves, temp_data)
}

# Create growth curves plot
growth_plot <- ggplot(growth_curves, aes(x = time, y = population, 
                                         color = species, linetype = model)) +
    geom_line(size = 1) +
    scale_y_log10(labels = scales::trans_format("log10", scales::math_format(10^.x))) +
    labs(title = "Bacterial Growth Patterns",
         subtitle = "Exponential vs Logistic Growth Models",
         x = "Time (hours)",
         y = "Population Size (log scale)",
         color = "Species",
         linetype = "Growth Model") +
    theme_minimal() +
    theme(legend.position = "bottom") +
    geom_hline(yintercept = K, linetype = "dashed", alpha = 0.5) +
    annotate("text", x = 20, y = K*1.5, label = "Carrying Capacity", size = 3)

print(growth_plot)

# Create doubling time comparison
doubling_plot <- ggplot(bacteria_data, aes(x = reorder(species, doubling_time_min), 
                                           y = doubling_time_min, 
                                           fill = environment)) +
    geom_col(alpha = 0.8) +
    scale_y_log10() +
    scale_fill_viridis_d() +
    labs(title = "Bacterial Doubling Times",
         x = "Species",
         y = "Doubling Time (minutes, log scale)",
         fill = "Environment") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    geom_text(aes(label = paste(doubling_time_min, "min")), 
              vjust = -0.2, size = 3)

print(doubling_plot)

cat("\n=== Prokaryotic Growth Analysis ===\n")
cat("Doubling times:\n")
for(i in 1:nrow(bacteria_data)) {
    cat(sprintf("  %s: %.0f minutes (%.2f hours)\n", 
                bacteria_data$species[i],
                bacteria_data$doubling_time_min[i],
                bacteria_data$doubling_time_min[i]/60))
}

# Calculate populations after 12 hours
cat("\nPopulation after 12 hours (exponential growth):\n")
for(i in 1:3) {
    final_pop <- exponential_growth(N0, bacteria_data$growth_rate[i], 12)
    cat(sprintf("  %s: %s cells\n", 
                bacteria_data$species[i],
                formatC(final_pop, format = "e", digits = 2)))
}

# Growth rate comparison
fastest <- bacteria_data[which.min(bacteria_data$doubling_time_min), ]
slowest <- bacteria_data[which.max(bacteria_data$doubling_time_min), ]
cat("\nGrowth rate extremes:\n")
cat("  Fastest:", fastest$species, "(", fastest$doubling_time_min, "min)\n")
cat("  Slowest:", slowest$species, "(", slowest$doubling_time_min, "min)\n")
cat("  Ratio:", round(slowest$doubling_time_min / fastest$doubling_time_min), "×\n")

# R: Eukaryotic Organelle Analysis
# Comparing organelle characteristics and functions

# Organelle characteristics
organelles <- data.frame(
    organelle = c("Nucleus", "Mitochondrion", "Chloroplast", "ER", "Golgi", 
                  "Lysosome", "Peroxisome", "Ribosome", "Vacuole"),
    size_um = c(10, 2, 8, 5, 3, 0.5, 0.8, 0.025, 15),
    membrane_layers = c(2, 2, 2, 1, 1, 1, 1, 0, 1),
    has_dna = c("Yes", "Yes", "Yes", "No", "No", "No", "No", "No", "No"),
    primary_function = c("Control", "Energy", "Photosynthesis", "Synthesis", 
                        "Processing", "Digestion", "Detox", "Translation", "Storage"),
    cell_type = c("All", "All", "Plant", "All", "All", "Animal", "All", "All", "Plant"),
    evolutionary_origin = c("Eukaryotic", "Endosymbiont", "Endosymbiont", "Eukaryotic",
                           "Eukaryotic", "Eukaryotic", "Endosymbiont", "Prokaryotic", "Eukaryotic"),
    relative_abundance = c(1, 100, 50, 1000, 50, 300, 100, 10000, 1)
)

# Size comparison plot
size_plot <- ggplot(organelles, aes(x = reorder(organelle, size_um), 
                                   y = size_um, 
                                   fill = evolutionary_origin)) +
    geom_col(alpha = 0.8, color = "black", size = 0.3) +
    scale_y_log10() +
    scale_fill_viridis_d() +
    coord_flip() +
    labs(title = "Organelle Sizes and Evolutionary Origins",
         x = "Organelle",
         y = "Size (μm, log scale)",
         fill = "Evolutionary Origin") +
    theme_minimal() +
    geom_text(aes(label = paste(size_um, "μm")), 
              hjust = -0.1, size = 3)

print(size_plot)

# Abundance comparison
abundance_plot <- ggplot(organelles[organelles$organelle != "ER", ], 
                        aes(x = reorder(organelle, relative_abundance), 
                           y = relative_abundance, 
                           fill = cell_type)) +
    geom_col(alpha = 0.8) +
    scale_y_log10(labels = scales::trans_format("log10", scales::math_format(10^.x))) +
    scale_fill_viridis_d() +
    coord_flip() +
    labs(title = "Relative Organelle Abundance in Cells",
         x = "Organelle",
         y = "Relative Number per Cell (log scale)",
         fill = "Found in") +
    theme_minimal()

print(abundance_plot)

# Analysis summary
cat("\n=== Eukaryotic Organelle Analysis ===\n")

# Endosymbiotic organelles
endosymbiont_count <- sum(organelles$evolutionary_origin == "Endosymbiont")
total_organelles <- nrow(organelles)
cat("Endosymbiotic organelles:", endosymbiont_count, "of", total_organelles, 
    "(", round(endosymbiont_count/total_organelles*100, 1), "%)\n")

# DNA-containing organelles
dna_organelles <- organelles[organelles$has_dna == "Yes", ]
cat("\nDNA-containing organelles:\n")
for(i in 1:nrow(dna_organelles)) {
    cat("  ", dna_organelles$organelle[i], ":", dna_organelles$primary_function[i], "\n")
}

# Size range
cat("\nSize range:", min(organelles$size_um), "to", max(organelles$size_um), "μm\n")
cat("Size ratio (largest/smallest):", round(max(organelles$size_um)/min(organelles$size_um)), "×\n")

# Membrane organization
cat("\nMembrane organization:\n")
membrane_summary <- table(organelles$membrane_layers)
cat("  No membrane:", membrane_summary["0"], "organelles\n")
cat("  Single membrane:", membrane_summary["1"], "organelles\n")
cat("  Double membrane:", membrane_summary["2"], "organelles\n")

# Most and least abundant
most_abundant <- organelles[which.max(organelles$relative_abundance), ]
least_abundant <- organelles[organelles$relative_abundance > 0 & 
                           organelles$relative_abundance == min(organelles$relative_abundance[organelles$relative_abundance > 0]), ]
cat("\nAbundance extremes:\n")
cat("  Most abundant:", most_abundant$organelle, "(", most_abundant$relative_abundance, "per cell)\n")
cat("  Least abundant:", least_abundant$organelle[1], "(", least_abundant$relative_abundance[1], "per cell)\n")

# R: Cellular Energy Analysis
# ATP production and cellular energy economics

# ATP yield from different metabolic pathways
atp_pathways <- data.frame(
    pathway = c("Glycolysis", "Citric Acid Cycle", "Electron Transport", "Total Aerobic", "Fermentation"),
    location = c("Cytoplasm", "Mitochondria", "Mitochondria", "Whole cell", "Cytoplasm"),
    atp_yield = c(2, 2, 28, 32, 2),
    oxygen_required = c(FALSE, TRUE, TRUE, TRUE, FALSE),
    efficiency_percent = c(6.25, 6.25, 87.5, 100, 6.25)
)

# Energy comparison data
energy_systems <- data.frame(
    system = c("Cellular respiration", "Photosynthesis", "Car engine", 
               "Steam engine", "LED light", "Incandescent bulb"),
    efficiency_percent = c(38, 3, 25, 10, 50, 5),
    type = c("Biological", "Biological", "Mechanical", "Mechanical", 
             "Electrical", "Electrical")
)

# ATP pathway visualization
pathway_plot <- ggplot(atp_pathways[1:4, ], aes(x = reorder(pathway, atp_yield), 
                                                y = atp_yield, 
                                                fill = oxygen_required)) +
    geom_col(alpha = 0.8, color = "black", size = 0.3) +
    coord_flip() +
    scale_fill_manual(values = c("FALSE" = "orange", "TRUE" = "steelblue"),
                     labels = c("Anaerobic", "Aerobic")) +
    labs(title = "ATP Production by Metabolic Pathway",
         subtitle = "Net ATP yield per glucose molecule",
         x = "Pathway",
         y = "Net ATP Yield",
         fill = "Oxygen Requirement") +
    theme_minimal() +
    geom_text(aes(label = paste(atp_yield, "ATP")), 
              hjust = -0.1, size = 4)

print(pathway_plot)

# Energy efficiency comparison
efficiency_plot <- ggplot(energy_systems, aes(x = reorder(system, efficiency_percent), 
                                              y = efficiency_percent, 
                                              fill = type)) +
    geom_col(alpha = 0.8, color = "black", size = 0.3) +
    coord_flip() +
    scale_fill_viridis_d() +
    labs(title = "Energy Conversion Efficiency Comparison",
         subtitle = "Biological vs Human-made Systems",
         x = "Energy System",
         y = "Efficiency (%)",
         fill = "System Type") +
    theme_minimal() +
    geom_text(aes(label = paste(efficiency_percent, "%")), 
              hjust = -0.1, size = 3)

print(efficiency_plot)

# Metabolic rates for different cell types
metabolic_rates <- data.frame(
    cell_type = c("Brain neuron", "Heart muscle", "Liver", "Kidney", 
                  "Skeletal muscle", "Fat", "Skin", "Bone"),
    atp_consumption_relative = c(100, 90, 75, 70, 60, 10, 15, 12),
    oxygen_consumption = c("Very High", "Very High", "High", "High", 
                          "Variable", "Low", "Low", "Low")
)

# Metabolic rate plot
metabolic_plot <- ggplot(metabolic_rates, aes(x = reorder(cell_type, atp_consumption_relative), 
                                              y = atp_consumption_relative, 
                                              fill = oxygen_consumption)) +
    geom_col(alpha = 0.8) +
    coord_flip() +
    scale_fill_viridis_d() +
    labs(title = "Metabolic Rates of Different Cell Types",
         x = "Cell Type",
         y = "Relative ATP Consumption",
         fill = "O₂ Consumption") +
    theme_minimal()

print(metabolic_plot)

# Analysis summary
cat("\n=== Cellular Energetics Analysis ===\n")

# ATP production efficiency
total_atp <- atp_pathways[atp_pathways$pathway == "Total Aerobic", "atp_yield"]
glucose_energy <- 2870  # kJ/mol
atp_energy <- 30.5      # kJ/mol

theoretical_efficiency <- (total_atp * atp_energy / glucose_energy) * 100
cat("ATP production from glucose:\n")
cat(sprintf("  Theoretical efficiency: %.1f%% (%d ATP × %.1f kJ/mol ÷ %d kJ/mol)\n", 
            theoretical_efficiency, total_atp, atp_energy, glucose_energy))

# Aerobic vs anaerobic comparison
aerobic_yield <- atp_pathways[atp_pathways$pathway == "Total Aerobic", "atp_yield"]
anaerobic_yield <- atp_pathways[atp_pathways$pathway == "Fermentation", "atp_yield"]
cat("\nAerobic vs Anaerobic:\n")
cat("  Aerobic yield:", aerobic_yield, "ATP per glucose\n")
cat("  Anaerobic yield:", anaerobic_yield, "ATP per glucose\n")
cat("  Aerobic advantage:", aerobic_yield/anaerobic_yield, "× more efficient\n")

# Biological vs artificial efficiency
bio_systems <- energy_systems[energy_systems$type == "Biological", ]
artificial_systems <- energy_systems[energy_systems$type != "Biological", ]

cat("\nEfficiency comparison:\n")
cat("  Best biological system:", bio_systems[which.max(bio_systems$efficiency_percent), "system"],
    "(", max(bio_systems$efficiency_percent), "%)\n")
cat("  Best artificial system:", artificial_systems[which.max(artificial_systems$efficiency_percent), "system"],
    "(", max(artificial_systems$efficiency_percent), "%)\n")

# Metabolic diversity
cat("\nCellular metabolic diversity:\n")
max_metabolic <- metabolic_rates[which.max(metabolic_rates$atp_consumption_relative), ]
min_metabolic <- metabolic_rates[which.min(metabolic_rates$atp_consumption_relative), ]
cat("  Most active:", max_metabolic$cell_type, "(", max_metabolic$atp_consumption_relative, "relative units)\n")
cat("  Least active:", min_metabolic$cell_type, "(", min_metabolic$atp_consumption_relative, "relative units)\n")
cat("  Metabolic range:", round(max_metabolic$atp_consumption_relative/min_metabolic$atp_consumption_relative, 1), "×\n")

# R: Cell Biology Applications Analysis
# Market trends and research impact assessment

# Biotechnology market data (simplified)
biotech_markets <- data.frame(
    application = c("Cell Therapy", "Drug Discovery", "Diagnostics", "Vaccines", 
                   "Tissue Engineering", "Gene Therapy", "Biomanufacturing"),
    market_size_billion_usd = c(15.2, 65.2, 45.1, 35.8, 8.9, 12.4, 28.7),
    growth_rate_percent = c(12.5, 8.7, 6.9, 9.2, 15.2, 18.4, 11.3),
    cell_based = c(TRUE, TRUE, FALSE, TRUE, TRUE, TRUE, TRUE),
    maturity = c("Emerging", "Mature", "Mature", "Mature", "Emerging", "Emerging", "Growing")
)

# Research publication trends (hypothetical data)
research_trends <- data.frame(
    year = 2015:2024,
    stem_cells = c(2500, 2800, 3200, 3800, 4200, 4800, 5200, 5800, 6200, 6800),
    cancer_biology = c(8000, 8500, 9200, 9800, 10500, 11200, 11800, 12500, 13100, 13800),
    gene_editing = c(800, 1200, 1800, 2800, 4200, 6000, 7200, 8500, 9800, 11000),
    synthetic_biology = c(600, 800, 1100, 1500, 2000, 2600, 3200, 3900, 4800, 5500)
)

# Create market analysis plot
market_plot <- ggplot(biotech_markets, aes(x = market_size_billion_usd, 
                                           y = growth_rate_percent, 
                                           color = maturity,
                                           size = market_size_billion_usd)) +
    geom_point(alpha = 0.8) +
    geom_text_repel(aes(label = application), size = 3) +
    scale_color_viridis_d() +
    labs(title = "Cell Biology Applications Market Analysis",
         subtitle = "Market Size vs Growth Rate (2024)",
         x = "Market Size (Billion USD)",
         y = "Annual Growth Rate (%)",
         color = "Market Maturity",
         size = "Market Size") +
    theme_minimal() +
    guides(size = "none") +
    geom_hline(yintercept = 10, linetype = "dashed", alpha = 0.5) +
    annotate("text", x = 60, y = 10.5, label = "High growth threshold", size = 3)

print(market_plot)

# Research trends visualization
# Reshape data for plotting
research_long <- data.frame(
    year = rep(research_trends$year, 4),
    publications = c(research_trends$stem_cells, research_trends$cancer_biology,
                    research_trends$gene_editing, research_trends$synthetic_biology),
    field = rep(c("Stem Cells", "Cancer Biology", "Gene Editing", "Synthetic Biology"), 
               each = length(research_trends$year))
)

research_plot <- ggplot(research_long, aes(x = year, y = publications, color = field)) +
    geom_line(size = 1.2) +
    geom_point(size = 2) +
    scale_color_viridis_d() +
    scale_x_continuous(breaks = seq(2015, 2024, 2)) +
    labs(title = "Research Publication Trends in Cell Biology",
         subtitle = "Annual Publications by Field (2015-2024)",
         x = "Year",
         y = "Number of Publications",
         color = "Research Field") +
    theme_minimal() +
    theme(legend.position = "bottom")

print(research_plot)

# Cell-based therapy timeline
therapy_milestones <- data.frame(
    year = c(1990, 1995, 2003, 2012, 2017, 2020, 2024),
    milestone = c("First gene therapy", "Stem cell discovery", "Human genome", 
                  "iPSC therapies", "CAR-T approval", "COVID vaccines", "CRISPR therapies"),
    impact = c(6, 8, 9, 8, 9, 10, 9),
    category = c("Gene", "Stem Cell", "Genomics", "Stem Cell", "Immune", "Vaccine", "Gene")
)

timeline_plot <- ggplot(therapy_milestones, aes(x = year, y = impact, color = category)) +
    geom_point(size = 4, alpha = 0.8) +
    geom_line(aes(group = 1), alpha = 0.4, size = 1) +
    geom_text_repel(aes(label = milestone), size = 3) +
    scale_x_continuous(breaks = seq(1990, 2025, 5)) +
    scale_y_continuous(breaks = 1:10) +
    scale_color_viridis_d() +
    labs(title = "Cell-Based Therapy Milestones",
         subtitle = "Major Breakthroughs in Clinical Applications (1990-2024)",
         x = "Year",
         y = "Clinical Impact (1-10)",
         color = "Therapy Type") +
    theme_minimal()

print(timeline_plot)

# Analysis summary
cat("\n=== Cell Biology Applications Analysis ===\n")

# Market analysis
total_market <- sum(biotech_markets$market_size_billion_usd)
cell_based_market <- sum(biotech_markets[biotech_markets$cell_based, "market_size_billion_usd"])
cell_based_percentage <- (cell_based_market / total_market) * 100

cat("Biotechnology market insights:\n")
cat("  Total market size:", total_market, "billion USD\n")
cat("  Cell-based applications:", round(cell_based_percentage, 1), "%\n")

fastest_growing <- biotech_markets[which.max(biotech_markets$growth_rate_percent), ]
largest_market <- biotech_markets[which.max(biotech_markets$market_size_billion_usd), ]
cat("  Fastest growing:", fastest_growing$application, 
    "(", fastest_growing$growth_rate_percent, "%/year)\n")
cat("  Largest market:", largest_market$application, 
    "(", largest_market$market_size_billion_usd, "billion USD)\n")

# Research trends analysis
cat("\nResearch publication growth (2015-2024):\n")
fields <- c("stem_cells", "cancer_biology", "gene_editing", "synthetic_biology")
field_names <- c("Stem Cells", "Cancer Biology", "Gene Editing", "Synthetic Biology")

for(i in 1:length(fields)) {
    initial <- research_trends[[fields[i]]][1]
    final <- research_trends[[fields[i]]][length(research_trends[[fields[i]]])]
    growth_factor <- final / initial
    cat(sprintf("  %s: %.1f× increase (%d to %d publications)\n", 
                field_names[i], growth_factor, initial, final))
}

# Therapy development timeline
cat("\nTherapy development timeline:\n")
timeline_span <- max(therapy_milestones$year) - min(therapy_milestones$year)
avg_impact <- mean(therapy_milestones$impact)
cat("  Development span:", timeline_span, "years\n")
cat("  Average clinical impact:", round(avg_impact, 1), "/10\n")

# Identify highest impact therapies
high_impact <- therapy_milestones[therapy_milestones$impact >= 9, ]
cat("  Major breakthroughs (impact ≥ 9):\n")
for(i in 1:nrow(high_impact)) {
    cat("    ", high_impact$year[i], ":", high_impact$milestone[i], "\n")
}