1. Descripción del proyecto

Se realizaron 6 simulaciones de dinámica molecular (10 ns c/u, GPU RTX 5080) de la proteína de fusión RSV-F (PDB: 5EA3) en GROMACS 2025.3:

  • 5 complejos proteína–ligando (AMBER99sb-ildn + GAFF2 / ACPYPE)
  • 1 simulación apo (proteína sola, control)

Parámetros: caja cúbica con borde de 1.2 nm, TIP3P, PME, 310 K (v-rescale), dt = 2 fs, 10 ns producción.

SYS_DIR <- "/home/svirology/DINÁMICAS_MULTISITIO_FINALES/systems_final_simple"

systems <- tribble(
  ~Label,                    ~Site,   ~Ligand,          ~Dir,
  "Site O — Rimegepant",    "Site O",  "Rimegepant",   file.path(SYS_DIR, "00_site_O_RIMEGEPANT_SULFATE_complex"),
  "Site II — Chlorhexidine","Site II", "Chlorhexidine", file.path(SYS_DIR, "02_site_II_CHLORHEXIDINE_complex"),
  "Site III — Rimegepant",  "Site III","Rimegepant",   file.path(SYS_DIR, "03_site_III_RIMEGEPANT_SULFATE_complex"),
  "Site IV — Atovaquone",   "Site IV", "Atovaquone",   file.path(SYS_DIR, "04_site_IV_ATOVAQUONE_complex"),
  "Site V — Vibegron",      "Site V",  "Vibegron",     file.path(SYS_DIR, "06_site_V_VIBEGRON_complex"),
  "Apo (control)",          "Apo",     "—",            file.path(SYS_DIR, "apo_protein")
)
kable(systems[, c("Label","Ligand","Dir")], col.names = c("Sistema","Ligando","Directorio")) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)
Sistema Ligando Directorio
Site O — Rimegepant Rimegepant /home/svirology/DINÁMICAS_MULTISITIO_FINALES/systems_final_simple/00_site_O_RIMEGEPANT_SULFATE_complex
Site II — Chlorhexidine Chlorhexidine /home/svirology/DINÁMICAS_MULTISITIO_FINALES/systems_final_simple/02_site_II_CHLORHEXIDINE_complex
Site III — Rimegepant Rimegepant /home/svirology/DINÁMICAS_MULTISITIO_FINALES/systems_final_simple/03_site_III_RIMEGEPANT_SULFATE_complex
Site IV — Atovaquone Atovaquone /home/svirology/DINÁMICAS_MULTISITIO_FINALES/systems_final_simple/04_site_IV_ATOVAQUONE_complex
Site V — Vibegron Vibegron /home/svirology/DINÁMICAS_MULTISITIO_FINALES/systems_final_simple/06_site_V_VIBEGRON_complex
Apo (control) — /home/svirology/DINÁMICAS_MULTISITIO_FINALES/systems_final_simple/apo_protein

2. Funciones auxiliares

# Lector XVG robusto (según skill md-analysis)
read_xvg <- function(filepath, col_names = c("X", "Y")) {
  if (!file.exists(filepath)) return(NULL)
  lines    <- readLines(filepath, warn = FALSE)
  data_ln  <- lines[!grepl("^[@#]", lines) & nchar(trimws(lines)) > 0]
  if (length(data_ln) == 0) return(NULL)
  df <- tryCatch(
    read.table(text = paste(data_ln, collapse = "\n"), fill = TRUE, comment.char = ""),
    error = function(e) NULL
  )
  if (is.null(df) || ncol(df) < 2) return(NULL)
  df <- df[, 1:min(length(col_names), ncol(df)), drop = FALSE]
  colnames(df) <- col_names[1:ncol(df)]
  df[is.finite(df[[1]]) & is.finite(df[[2]]), ]
}

# Adjunta metadatos y convierte ps→ns si max(X) > 500
attach_meta <- function(df, label, time_col = "X") {
  if (is.null(df) || nrow(df) == 0) return(NULL)
  x_vals <- df[[time_col]]
  if (max(x_vals, na.rm = TRUE) > 500) x_vals <- x_vals / 1000
  df$Time_ns <- x_vals
  df$Label   <- label
  df
}

# Combina lista de data.frames ignorando NULLs
safe_bind <- function(lst) {
  items <- Filter(Negate(is.null), lst)
  if (length(items) == 0) return(data.frame())
  do.call(rbind, items)
}

3. Carga de datos

all_data <- list()

for (i in seq_len(nrow(systems))) {
  s     <- systems[i, ]
  d     <- s$Dir
  label <- s$Label
  entry <- list(label = label, site = s$Site, ligand = s$Ligand)

  # RMSD backbone
  df <- read_xvg(file.path(d, "rmsd_backbone.xvg"))
  if (!is.null(df)) {
    df$Y <- df$Y * 10   # nm → Å
    entry$rmsd <- attach_meta(df, label)
  }

  # RMSD ligando (solo complejos)
  df <- read_xvg(file.path(d, "rmsd_ligand.xvg"))
  if (!is.null(df)) {
    df$Y <- df$Y * 10
    entry$rmsd_lig <- attach_meta(df, label)
  }

  # RMSF
  df <- read_xvg(file.path(d, "rmsf_residue.xvg"), c("Residue","RMSF"))
  if (!is.null(df)) {
    df$RMSF  <- df$RMSF * 10   # nm → Å
    df$Label <- label
    entry$rmsf <- df
  }

  # Energía potencial (en ps, attach_meta convierte automáticamente)
  df <- read_xvg(file.path(d, "energy.xvg"), c("X","Potential"))
  if (!is.null(df)) {
    entry$energy <- attach_meta(df, label)
  }

  # SASA
  df <- read_xvg(file.path(d, "sasa.xvg"), c("X","SASA"))
  if (!is.null(df)) {
    entry$sasa <- attach_meta(df, label)
  }

  # Radio de giro
  df <- read_xvg(file.path(d, "gyrate.xvg"), c("X","Rg"))
  if (!is.null(df)) {
    df$Rg <- df$Rg * 10   # nm → Å
    entry$gyrate <- attach_meta(df, label)
  }

  # H-bonds desde caché numpy (via Python → CSV temporal)
  hb_cache <- file.path(d, "hbonds_mda.npy")
  if (file.exists(hb_cache) && s$Ligand != "—") {
    csv_tmp <- file.path(d, "hbonds_mda.csv")
    if (!file.exists(csv_tmp)) {
      system(paste0(
        "python3 -c \"import numpy as np; arr=np.load('", hb_cache, "'); ",
        "import pandas as pd; ",
        "pd.DataFrame({'frame':range(len(arr)),'hbonds':arr}).to_csv('", csv_tmp, "',index=False)\""
      ))
    }
    if (file.exists(csv_tmp)) {
      hb <- read.csv(csv_tmp)
      hb$Time_ns <- seq(0, 10, length.out = nrow(hb))
      hb$Label   <- label
      entry$hbonds <- hb
    }
  }

  all_data[[label]] <- entry
}

cat("Sistemas cargados:", length(all_data), "\n")
## Sistemas cargados: 6
cat("Métricas disponibles por sistema:\n")
## Métricas disponibles por sistema:
for (nm in names(all_data)) {
  metricas <- names(all_data[[nm]])[!names(all_data[[nm]]) %in% c("label","site","ligand")]
  cat(sprintf("  %-30s %s\n", nm, paste(metricas, collapse=", ")))
}
##   Site O — Rimegepant          rmsd, rmsd_lig, rmsf, energy, sasa, gyrate, hbonds
##   Site II — Chlorhexidine      rmsd, rmsd_lig, rmsf, energy, sasa, gyrate, hbonds
##   Site III — Rimegepant        rmsd, rmsd_lig, rmsf, energy, sasa, gyrate, hbonds
##   Site IV — Atovaquone         rmsd, rmsd_lig, rmsf, energy, sasa, gyrate, hbonds
##   Site V — Vibegron            rmsd, rmsd_lig, rmsf, energy, sasa, gyrate, hbonds
##   Apo (control)                  rmsd, rmsf, energy, sasa, gyrate

4. RMSD Backbone (estabilidad estructural)

rmsd_df <- safe_bind(lapply(all_data, function(x) x$rmsd))

p_bb <- ggplot(rmsd_df, aes(x = Time_ns, y = Y, color = Label, linetype = Label)) +
  geom_line(alpha = 0.85, size = 0.7) +
  scale_color_manual(values = site_colors) +
  scale_linetype_manual(values = c(
    "Site O — Rimegepant"    = "solid",
    "Site II — Chlorhexidine"= "solid",
    "Site III — Rimegepant"  = "solid",
    "Site IV — Atovaquone"   = "solid",
    "Site V — Vibegron"      = "solid",
    "Apo (control)"          = "dashed"
  )) +
  labs(title    = "RMSD Backbone Cα — RSV-F (5EA3)",
       subtitle = "Referencia: estructura inicial minimizada",
       x = "Tiempo (ns)", y = "RMSD (Å)",
       color = "Sistema", linetype = "Sistema") +
  theme_bw(base_size = 12) +
  theme(legend.position = "right") +
  coord_cartesian(xlim = c(0, 10))

print(p_bb)

rmsd_stats <- rmsd_df %>%
  group_by(Label) %>%
  summarise(
    Media_Å  = round(mean(Y, na.rm = TRUE), 2),
    SD_Å     = round(sd(Y, na.rm = TRUE), 2),
    Min_Å    = round(min(Y, na.rm = TRUE), 2),
    Max_Å    = round(max(Y, na.rm = TRUE), 2),
    .groups = "drop"
  )

rmsd_stats %>%
  kable(caption = "Estadísticas RMSD Backbone (Å)") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE) %>%
  row_spec(which(rmsd_stats$Media_Å == min(rmsd_stats$Media_Å)), bold = TRUE, color = "white", background = "#2A9D8F")
Estadísticas RMSD Backbone (Å)
Label Media_Å SD_Å Min_Å Max_Å
Apo (control) 2.86 0.59 0.01 4.77
Site II — Chlorhexidine 3.60 0.58 0.01 4.74
Site III — Rimegepant 2.01 0.30 0.01 2.90
Site IV — Atovaquone 2.77 0.49 0.00 4.04
Site O — Rimegepant 2.46 0.43 0.01 3.97
Site V — Vibegron 2.50 0.66 0.01 4.26

5. RMSD Ligando (estabilidad de unión)

rmsd_lig_df <- safe_bind(lapply(all_data, function(x) x$rmsd_lig))

# Panel A: todos los sitios
p_lig_all <- ggplot(rmsd_lig_df, aes(x = Time_ns, y = Y, color = Label)) +
  geom_line(alpha = 0.85, size = 0.7) +
  scale_color_manual(values = site_colors) +
  labs(title = "A) RMSD Ligando — todos los sitios",
       x = "Tiempo (ns)", y = "RMSD (Å)", color = "Sistema") +
  theme_bw(base_size = 12) +
  coord_cartesian(xlim = c(0, 10))

# Panel B: zoom sitios estables (sin Site O)
rmsd_lig_stable <- rmsd_lig_df %>%
  filter(!grepl("Site O", Label))

p_lig_stable <- ggplot(rmsd_lig_stable, aes(x = Time_ns, y = Y, color = Label)) +
  geom_line(alpha = 0.85, size = 0.7) +
  scale_color_manual(values = site_colors) +
  labs(title = "B) RMSD Ligando — Sitios II–V (zoom)",
       x = "Tiempo (ns)", y = "RMSD (Å)", color = "Sistema") +
  theme_bw(base_size = 12) +
  coord_cartesian(xlim = c(0, 10))

p_lig_all / p_lig_stable

lig_stats <- rmsd_lig_df %>%
  group_by(Label) %>%
  summarise(
    Media_Å  = round(mean(Y, na.rm = TRUE), 2),
    SD_Å     = round(sd(Y, na.rm = TRUE), 2),
    Max_Å    = round(max(Y, na.rm = TRUE), 2),
    Estabilidad = case_when(
      mean(Y, na.rm = TRUE) < 3  ~ "Muy estable ✅",
      mean(Y, na.rm = TRUE) < 8  ~ "Estable ✅",
      mean(Y, na.rm = TRUE) < 20 ~ "Moderado ⚠️",
      TRUE                       ~ "Escapó ❌"
    ),
    .groups = "drop"
  )

lig_stats %>%
  kable(caption = "Estadísticas RMSD Ligando (Å)") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE) %>%
  row_spec(which(lig_stats$Media_Å == min(lig_stats$Media_Å)), bold = TRUE,
           color = "white", background = "#2A9D8F") %>%
  row_spec(which(lig_stats$Media_Å == max(lig_stats$Media_Å)), bold = TRUE,
           color = "white", background = "#E63946")
Estadísticas RMSD Ligando (Å)
Label Media_Å SD_Å Max_Å Estabilidad
Site II — Chlorhexidine 3.20 0.88 6.62 Estable ✅
Site III — Rimegepant 1.75 0.75 7.78 Muy estable ✅
Site IV — Atovaquone 7.26 4.32 16.82 Estable ✅
Site O — Rimegepant 69.86 37.42 122.71 Escapó ❌
Site V — Vibegron 6.89 2.46 12.48 Estable ✅

6. RMSF por residuo (flexibilidad)

rmsf_df <- safe_bind(lapply(all_data, function(x) x$rmsf))

# Panel por sistema
ggplot(rmsf_df, aes(x = Residue, y = RMSF, color = Label, fill = Label)) +
  geom_line(size = 0.7, alpha = 0.9) +
  geom_area(alpha = 0.12, position = "identity") +
  scale_color_manual(values = site_colors) +
  scale_fill_manual(values = site_colors) +
  facet_wrap(~ Label, ncol = 2, scales = "free_y") +
  labs(title    = "RMSF por Residuo — Flexibilidad del Backbone",
       subtitle = "Picos indican regiones de alta movilidad (loops, extremos)",
       x = "Número de residuo", y = "RMSF (Å)") +
  theme_bw(base_size = 11) +
  theme(legend.position = "none",
        strip.background = element_rect(fill = "#264653"),
        strip.text = element_text(color = "white", face = "bold"))

# Overlay comparativo
ggplot(rmsf_df, aes(x = Residue, y = RMSF, color = Label)) +
  geom_line(size = 0.6, alpha = 0.8) +
  scale_color_manual(values = site_colors) +
  labs(title = "RMSF Overlay — comparación de todos los sistemas",
       x = "Residuo", y = "RMSF (Å)", color = "Sistema") +
  theme_bw(base_size = 12)

7. Radio de giro (compacidad estructural)

gyrate_df <- safe_bind(lapply(all_data, function(x) x$gyrate))

ggplot(gyrate_df, aes(x = Time_ns, y = Rg, color = Label, linetype = Label)) +
  geom_line(alpha = 0.85, size = 0.7) +
  scale_color_manual(values = site_colors) +
  scale_linetype_manual(values = c(
    "Site O — Rimegepant"    = "solid",
    "Site II — Chlorhexidine"= "solid",
    "Site III — Rimegepant"  = "solid",
    "Site IV — Atovaquone"   = "solid",
    "Site V — Vibegron"      = "solid",
    "Apo (control)"          = "dashed"
  )) +
  labs(title    = "Radio de Giro (Rg) — Compacidad estructural",
       subtitle = "Rg estable indica que la proteína no se despliega",
       x = "Tiempo (ns)", y = "Rg (Å)",
       color = "Sistema", linetype = "Sistema") +
  theme_bw(base_size = 12) +
  coord_cartesian(xlim = c(0, 10))

gyrate_df %>%
  group_by(Label) %>%
  summarise(
    Rg_medio_Å = round(mean(Rg, na.rm = TRUE), 2),
    Rg_SD_Å    = round(sd(Rg, na.rm = TRUE), 2),
    Rg_min_Å   = round(min(Rg, na.rm = TRUE), 2),
    Rg_max_Å   = round(max(Rg, na.rm = TRUE), 2),
    .groups = "drop"
  ) %>%
  kable(caption = "Estadísticas Radio de Giro (Å)") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)
Estadísticas Radio de Giro (Å)
Label Rg_medio_Å Rg_SD_Å Rg_min_Å Rg_max_Å
Apo (control) 29.05 0.39 27.96 30.10
Site II — Chlorhexidine 29.06 0.39 27.67 29.97
Site III — Rimegepant 28.86 0.22 28.17 29.52
Site IV — Atovaquone 28.14 0.44 27.13 29.52
Site O — Rimegepant 28.47 0.42 27.34 29.49
Site V — Vibegron 28.83 0.35 27.96 29.97

8. Energía potencial

energy_df <- safe_bind(lapply(all_data, function(x) x$energy))

ggplot(energy_df, aes(x = Time_ns, y = Potential / 1000, color = Label, linetype = Label)) +
  geom_line(alpha = 0.8, size = 0.6) +
  scale_color_manual(values = site_colors) +
  scale_linetype_manual(values = c(
    "Site O — Rimegepant"    = "solid",
    "Site II — Chlorhexidine"= "solid",
    "Site III — Rimegepant"  = "solid",
    "Site IV — Atovaquone"   = "solid",
    "Site V — Vibegron"      = "solid",
    "Apo (control)"          = "dashed"
  )) +
  labs(title    = "Energía Potencial del sistema",
       subtitle = "Valores estables indican simulaciones bien equilibradas",
       x = "Tiempo (ns)", y = "Energía potencial (kJ/mol × 10³)",
       color = "Sistema", linetype = "Sistema") +
  theme_bw(base_size = 12) +
  coord_cartesian(xlim = c(0, 10))

energy_df %>%
  group_by(Label) %>%
  summarise(
    E_media   = round(mean(Potential, na.rm = TRUE)),
    E_SD      = round(sd(Potential, na.rm = TRUE)),
    E_drift   = round(last(Potential) - first(Potential)),
    .groups = "drop"
  ) %>%
  kable(col.names = c("Sistema","E media (kJ/mol)","SD (kJ/mol)","Drift total (kJ/mol)"),
        caption = "Energía Potencial — estadísticas de simulación") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE)
Energía Potencial — estadísticas de simulación
Sistema E media (kJ/mol) SD (kJ/mol) Drift total (kJ/mol)
Apo (control) -2521840 1642 3482
Site II — Chlorhexidine -2522524 1634 -580
Site III — Rimegepant -2521807 1658 -357
Site IV — Atovaquone -2521712 1591 1117
Site O — Rimegepant -2590529 1623 -1365
Site V — Vibegron -2522930 1673 -3052

9. SASA — Superficie accesible al solvente

sasa_df <- safe_bind(lapply(all_data, function(x) x$sasa))

ggplot(sasa_df, aes(x = Time_ns, y = SASA, color = Label, linetype = Label)) +
  geom_line(alpha = 0.85, size = 0.7) +
  scale_color_manual(values = site_colors) +
  scale_linetype_manual(values = c(
    "Site O — Rimegepant"    = "solid",
    "Site II — Chlorhexidine"= "solid",
    "Site III — Rimegepant"  = "solid",
    "Site IV — Atovaquone"   = "solid",
    "Site V — Vibegron"      = "solid",
    "Apo (control)"          = "dashed"
  )) +
  labs(title    = "SASA — Superficie Accesible al Solvente",
       subtitle = "Cambios en SASA reflejan apertura/cierre de cavidades",
       x = "Tiempo (ns)", y = "SASA (nm²)",
       color = "Sistema", linetype = "Sistema") +
  theme_bw(base_size = 12) +
  coord_cartesian(xlim = c(0, 10))

10. Puentes de hidrógeno proteína–ligando

hbond_list <- lapply(all_data, function(x) x$hbonds)
hbond_df   <- safe_bind(Filter(Negate(is.null), hbond_list))

if (nrow(hbond_df) > 0) {
  # Suavizado: media móvil de 50 frames
  smooth_hb <- function(x, w = 50) {
    if (length(x) < w) return(x)
    stats::filter(x, rep(1/w, w), sides = 2)
  }

  hbond_smooth <- hbond_df %>%
    group_by(Label) %>%
    mutate(hb_smooth = as.numeric(smooth_hb(hbonds, 50))) %>%
    ungroup()

  ggplot(hbond_smooth, aes(x = Time_ns, y = hb_smooth, color = Label)) +
    geom_line(alpha = 0.9, size = 0.9, na.rm = TRUE) +
    scale_color_manual(values = site_colors) +
    labs(title    = "Puentes de Hidrógeno Proteína–Ligando (MDAnalysis)",
         subtitle = "Criterio geométrico: d(D-A) ≤ 3.5 Å, ángulo D-H-A ≥ 150°",
         x = "Tiempo (ns)", y = "N° H-bonds (media móvil 50 frames)",
         color = "Sistema") +
    theme_bw(base_size = 12) +
    coord_cartesian(xlim = c(0, 10), ylim = c(0, NA))
}

if (nrow(hbond_df) > 0) {
  hbond_df %>%
    group_by(Label) %>%
    summarise(
      HBonds_medio = round(mean(hbonds, na.rm = TRUE), 3),
      HBonds_max   = max(hbonds, na.rm = TRUE),
      Frames_con_HB = sum(hbonds > 0),
      Pct_frames   = round(100 * mean(hbonds > 0), 1),
      .groups = "drop"
    ) %>%
    arrange(desc(HBonds_medio)) %>%
    kable(col.names = c("Sistema","H-bonds medio","H-bonds máx","Frames con HB","% frames"),
          caption = "Estadísticas de Puentes de Hidrógeno") %>%
    kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE) %>%
    row_spec(1, bold = TRUE, color = "white", background = "#F4A261")
}
Estadísticas de Puentes de Hidrógeno
Sistema H-bonds medio H-bonds máx Frames con HB % frames
Site V — Vibegron 1.357 3 947 94.6
Site III — Rimegepant 0.864 3 837 83.6
Site IV — Atovaquone 0.261 2 232 23.2
Site O — Rimegepant 0.018 1 18 1.8
Site II — Chlorhexidine 0.000 0 0 0.0

11. Mapa de contactos proteína–ligando en el tiempo

library(dplyr); library(tidyr); library(ggplot2); library(patchwork)

CONTACTS_DIR <- "/home/svirology/DINÁMICAS_MULTISITIO_FINALES/analisis_resultados"

contact_files <- list(
  "Site O — Rimegepant"    = "contacts_time_00_site_O_RIMEGEPANT_SULFATE_complex.csv",
  "Site II — Chlorhexidine"= "contacts_time_02_site_II_CHLORHEXIDINE_complex.csv",
  "Site III — Rimegepant"  = "contacts_time_03_site_III_RIMEGEPANT_SULFATE_complex.csv",
  "Site IV — Atovaquone"   = "contacts_time_04_site_IV_ATOVAQUONE_complex.csv",
  "Site V — Vibegron"      = "contacts_time_06_site_V_VIBEGRON_complex.csv"
)

# Leer y combinar
contacts_all <- lapply(names(contact_files), function(sys_label) {
  path <- file.path(CONTACTS_DIR, contact_files[[sys_label]])
  if (!file.exists(path)) return(NULL)
  df <- read.csv(path, stringsAsFactors = FALSE)
  df$System <- sys_label
  df
}) %>% Filter(Negate(is.null), .) %>% bind_rows()

# Etiqueta residuo limpia
contacts_all <- contacts_all %>%
  mutate(res_label = paste0(resname, resid))

cat("Sistemas cargados:", n_distinct(contacts_all$System), "\n")
## Sistemas cargados: 5
cat("Total registros:", nrow(contacts_all), "\n")
## Total registros: 2476
# Función: heatmap + curva de contactos totales para un sistema
plot_contact_heatmap <- function(sys_label, top_n = 20) {
  df <- contacts_all %>% filter(System == sys_label)
  if (nrow(df) == 0) return(NULL)

  # Seleccionar los top_n residuos más contactados (total acumulado)
  top_res <- df %>%
    group_by(res_label) %>%
    summarise(total = sum(n_frames_contact), .groups = "drop") %>%
    slice_max(total, n = top_n) %>%
    pull(res_label)

  df_top <- df %>% filter(res_label %in% top_res)

  # Ordenar residuos por posición en secuencia (resid numérico)
  res_order <- df_top %>%
    distinct(res_label, resid) %>%
    arrange(resid) %>%
    pull(res_label)

  df_top <- df_top %>%
    mutate(res_label = factor(res_label, levels = res_order))

  # ── Heatmap ──────────────────────────────────────────────────────────────
  p_heat <- ggplot(df_top,
                   aes(x = time_ns, y = res_label, fill = contact_fraction)) +
    geom_tile(color = NA) +
    scale_fill_gradientn(
      colours = c("#f7fbff","#c6dbef","#6baed6","#2171b5","#08306b"),
      limits  = c(0, 1),
      labels  = scales::percent,
      name    = "Frames\ncon contacto"
    ) +
    scale_x_continuous(breaks = seq(0, 10, 2), expand = c(0, 0)) +
    labs(
      title    = sys_label,
      subtitle = paste0("Top ", top_n,
                        " residuos · Corte 4.5 Å · Ventanas de 0.2 ns"),
      x = "Tiempo (ns)",
      y = NULL
    ) +
    theme_bw(base_size = 11) +
    theme(
      axis.text.y     = element_text(size = 7.5, family = "mono"),
      legend.position = "right",
      panel.grid      = element_blank(),
      plot.title      = element_text(face = "bold", size = 12),
      plot.subtitle   = element_text(size = 9, color = "grey40")
    )

  # ── Curva de contactos totales por ventana ───────────────────────────────
  curve_df <- df %>%
    group_by(time_ns, bin) %>%
    summarise(
      n_res_contacted   = n_distinct(res_label),
      mean_frac_contact = mean(contact_fraction),
      .groups = "drop"
    )

  p_curve <- ggplot(curve_df, aes(x = time_ns)) +
    geom_ribbon(aes(ymin = 0, ymax = n_res_contacted),
                fill = "#2171b5", alpha = 0.15) +
    geom_line(aes(y = n_res_contacted, color = "Residuos en contacto"),
              size = 0.9) +
    geom_line(aes(y = mean_frac_contact * max(n_res_contacted),
                  color = "Fracción media contacto (escala)"),
              size = 0.7, linetype = "dashed") +
    scale_color_manual(
      values = c("Residuos en contacto"                = "#08306b",
                 "Fracción media contacto (escala)" = "#D62728"),
      name = NULL
    ) +
    scale_x_continuous(breaks = seq(0, 10, 2), expand = c(0, 0)) +
    scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
    labs(x = "Tiempo (ns)", y = "N° residuos") +
    theme_bw(base_size = 10) +
    theme(legend.position  = "bottom",
          legend.text      = element_text(size = 8),
          panel.grid.minor = element_blank())

  p_heat / p_curve + plot_layout(heights = c(3, 1))
}
plot_contact_heatmap("Site O — Rimegepant")

plot_contact_heatmap("Site II — Chlorhexidine")

plot_contact_heatmap("Site III — Rimegepant")

plot_contact_heatmap("Site IV — Atovaquone")

plot_contact_heatmap("Site V — Vibegron")

# Tabla: residuos más contactados por sistema
contacts_top_table <- contacts_all %>%
  group_by(System, res_label) %>%
  summarise(
    total_frames  = sum(n_frames_contact),
    mean_fraction = round(mean(contact_fraction), 3),
    n_windows     = n_distinct(bin),
    windows_present = sum(contact_fraction > 0),
    pct_time      = round(100 * windows_present / n_windows, 1),
    .groups = "drop"
  ) %>%
  group_by(System) %>%
  slice_max(mean_fraction, n = 5) %>%
  arrange(System, desc(mean_fraction))

contacts_top_table %>%
  select(System, res_label, mean_fraction, pct_time) %>%
  rename(Sistema = System,
         Residuo = res_label,
         `Fracción media` = mean_fraction,
         `% tiempo en contacto` = pct_time) %>%
  kable(caption = "Top 5 residuos más contactados por sistema (4.5 Å, 10 ns)") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE) %>%
  collapse_rows(columns = 1, valign = "top")
Top 5 residuos más contactados por sistema (4.5 Å, 10 ns)
Sistema Residuo Fracción media % tiempo en contacto
Site II — Chlorhexidine ASN181 1.000 100
LYS184 1.000 100
LYS185 1.000 100
ASP182 0.998 100
SER188 0.997 100
Site III — Rimegepant ALA268 1.000 100
ARG252 0.999 100
GLY253 0.999 100
LEU229 0.997 100
HIS230 0.996 100
Site IV — Atovaquone SER338 0.950 100
LYS340 0.923 100
LEU369 0.909 100
ILE344 0.892 100
SER364 0.876 100
Site O — Rimegepant LYS39 0.795 100
ILE38 0.743 100
ASN37 0.583 100
LYS40 0.435 100
VAL41 0.431 100
Site V — Vibegron ASN123 1.000 100
TYR124 1.000 100
LEU98 0.994 100
LYS94 0.951 100
ASN91 0.938 100

12. Poses del complejo proteína–ligando a t = 0, 5 y 10 ns

Estructuras extraídas directamente de la trayectoria de MD. Azul/verde/púrpura: proteína (cartoon). Naranja brillante: residuos en contacto ≤ 4.5 Å. Color sólido: ligando (sticks + esferas). Sin solvente ni iones.

POSES_DIR <- "/home/svirology/DINÁMICAS_MULTISITIO_FINALES/analisis_resultados/poses_timepoints"

pose_systems <- list(
  list(site="Site_O",   ligand="Rimegepant",    label="Site Ø — Rimegepant"),
  list(site="Site_II",  ligand="Chlorhexidine", label="Site II — Chlorhexidine"),
  list(site="Site_III", ligand="Rimegepant",    label="Site III — Rimegepant"),
  list(site="Site_IV",  ligand="Atovaquone",    label="Site IV — Atovaquone"),
  list(site="Site_V",   ligand="Vibegron",      label="Site V — Vibegron")
)
library(knitr)

for (sys in pose_systems) {
  imgs <- lapply(c(0, 5, 10), function(t) {
    file.path(POSES_DIR,
              paste0(sys$site, "_", sys$ligand, "_t", t, "ns.png"))
  })
  exist <- sapply(imgs, file.exists)

  if (!all(exist)) {
    cat(paste0("\n> ⚠ Imágenes no encontradas para ", sys$label, "\n\n"))
    next
  }

  cat(paste0("\n### ", sys$label, "\n\n"))

  # Componer panel 3 columnas con knitr::include_graphics
  cat('<div style="display:flex; gap:6px; align-items:flex-start; margin-bottom:18px;">\n')
  for (i in seq_along(imgs)) {
    t_label <- c("t = 0 ns", "t = 5 ns", "t = 10 ns")[i]
    cat(paste0(
      '<div style="flex:1; text-align:center;">',
      '<p style="font-size:12px; font-weight:bold; margin:0 0 4px 0; color:#444;">',
      t_label, '</p>',
      '<img src="', imgs[[i]], '" style="width:100%; border:1px solid #ddd;',
      ' border-radius:4px;" />',
      '</div>\n'
    ))
  }
  cat('</div>\n\n')
}

Site Ø — Rimegepant

t = 0 ns

t = 5 ns

t = 10 ns

Site II — Chlorhexidine

t = 0 ns

t = 5 ns

t = 10 ns

Site III — Rimegepant

t = 0 ns

t = 5 ns

t = 10 ns

Site IV — Atovaquone

t = 0 ns

t = 5 ns

t = 10 ns

Site V — Vibegron

t = 0 ns

t = 5 ns

t = 10 ns

13. Panel integrado de resultados

# Preparar subsets
rmsd_bb_5  <- rmsd_df %>% filter(Label != "Apo (control)")
rmsd_bb_apo <- rmsd_df %>% filter(Label == "Apo (control)")

p1 <- ggplot(rmsd_df, aes(x = Time_ns, y = Y, color = Label, linetype = Label)) +
  geom_line(alpha = 0.85, size = 0.65) +
  scale_color_manual(values = site_colors) +
  scale_linetype_manual(values = setNames(
    c("solid","solid","solid","solid","solid","dashed"), names(site_colors))) +
  labs(title = "A) RMSD Backbone Cα", x = "Tiempo (ns)", y = "RMSD (Å)") +
  theme_bw(base_size = 11) + coord_cartesian(xlim = c(0,10)) +
  theme(legend.position = "none")

p2 <- ggplot(rmsd_lig_df %>% filter(!grepl("Site O", Label)),
             aes(x = Time_ns, y = Y, color = Label)) +
  geom_line(alpha = 0.85, size = 0.65) +
  geom_line(data = rmsd_lig_df %>% filter(grepl("Site O", Label)),
            aes(x = Time_ns, y = Y, color = Label), linetype = "dashed", size = 0.65) +
  scale_color_manual(values = site_colors) +
  labs(title = "B) RMSD Ligando", x = "Tiempo (ns)", y = "RMSD (Å)") +
  theme_bw(base_size = 11) + coord_cartesian(xlim = c(0,10)) +
  theme(legend.position = "none")

p3 <- ggplot(gyrate_df, aes(x = Time_ns, y = Rg, color = Label, linetype = Label)) +
  geom_line(alpha = 0.85, size = 0.65) +
  scale_color_manual(values = site_colors) +
  scale_linetype_manual(values = setNames(
    c("solid","solid","solid","solid","solid","dashed"), names(site_colors))) +
  labs(title = "C) Radio de Giro", x = "Tiempo (ns)", y = "Rg (Å)") +
  theme_bw(base_size = 11) + coord_cartesian(xlim = c(0,10)) +
  theme(legend.position = "none")

p4 <- ggplot(sasa_df, aes(x = Time_ns, y = SASA, color = Label, linetype = Label)) +
  geom_line(alpha = 0.85, size = 0.65) +
  scale_color_manual(values = site_colors) +
  scale_linetype_manual(values = setNames(
    c("solid","solid","solid","solid","solid","dashed"), names(site_colors))) +
  labs(title = "D) SASA", x = "Tiempo (ns)", y = "SASA (nm²)") +
  theme_bw(base_size = 11) + coord_cartesian(xlim = c(0,10)) +
  theme(legend.position = "bottom",
        legend.text = element_text(size = 8))

(p1 | p2) / (p3 | p4) +
  plot_annotation(
    title    = "Dinámicas Moleculares RSV-F (5EA3) — 5 sitios de unión + control apo (10 ns, GPU)",
    subtitle = "AMBER99sb-ildn + GAFF2 | TIP3P | 310 K | dt=2 fs",
    theme    = theme(plot.title = element_text(size = 14, face = "bold"),
                     plot.subtitle = element_text(size = 10))
  )

14. Tabla resumen comparativa

# Construir tabla resumen completa
resumen <- systems %>%
  rowwise() %>%
  mutate(
    RMSD_prot_A = {
      d <- all_data[[Label]]
      if (!is.null(d$rmsd)) round(mean(d$rmsd$Y, na.rm=TRUE), 2) else NA
    },
    RMSD_lig_A = {
      d <- all_data[[Label]]
      if (!is.null(d$rmsd_lig)) round(mean(d$rmsd_lig$Y, na.rm=TRUE), 2) else NA
    },
    Rg_A = {
      d <- all_data[[Label]]
      if (!is.null(d$gyrate)) round(mean(d$gyrate$Rg, na.rm=TRUE), 2) else NA
    },
    SASA_nm2 = {
      d <- all_data[[Label]]
      if (!is.null(d$sasa)) round(mean(d$sasa$SASA, na.rm=TRUE), 2) else NA
    },
    HBonds = {
      d <- all_data[[Label]]
      if (!is.null(d$hbonds)) round(mean(d$hbonds$hbonds, na.rm=TRUE), 3) else NA
    },
    Estabilidad = {
      rl <- all_data[[Label]]$rmsd_lig
      if (is.null(rl)) "—"
      else {
        m <- mean(rl$Y, na.rm=TRUE)
        if (m < 3) "Muy estable ✅" else if (m < 8) "Estable ✅" else if (m < 20) "Moderado ⚠️" else "Escapó ❌"
      }
    }
  ) %>%
  ungroup() %>%
  select(Label, Ligand, RMSD_prot_A, RMSD_lig_A, Rg_A, SASA_nm2, HBonds, Estabilidad)

resumen %>%
  kable(
    col.names = c("Sistema","Ligando","RMSD prot (Å)","RMSD lig (Å)",
                  "Rg (Å)","SASA (nm²)","H-bonds","Estabilidad"),
    caption = "Resumen comparativo — Dinámica Molecular RSV-F 10 ns"
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = TRUE) %>%
  row_spec(which(resumen$Label == "Site III — Rimegepant"), bold = TRUE,
           background = "#fff3cd") %>%
  row_spec(which(resumen$Label == "Apo (control)"), italic = TRUE,
           background = "#f0f0f0")
Resumen comparativo — Dinámica Molecular RSV-F 10 ns
Sistema Ligando RMSD prot (Å) RMSD lig (Å) Rg (Å) SASA (nm²) H-bonds Estabilidad
Site O — Rimegepant Rimegepant 2.46 69.86 28.47 214.80 0.018 Escapó ❌
Site II — Chlorhexidine Chlorhexidine 3.60 3.20 29.06 219.05 0.000 Estable ✅
Site III — Rimegepant Rimegepant 2.01 1.75 28.86 214.20 0.864 Muy estable ✅
Site IV — Atovaquone Atovaquone 2.77 7.26 28.14 216.59 0.261 Estable ✅
Site V — Vibegron Vibegron 2.50 6.89 28.83 217.56 1.357 Estable ✅
Apo (control) — 2.86 NA 29.05 216.36 NA —

15. Ranking de candidatos

# Score compuesto (mayor es mejor):
# RMSD ligando inverso + H-bonds
candidatos <- resumen %>%
  filter(Ligand != "—", !is.na(RMSD_lig_A)) %>%
  mutate(
    Score_RMSD  = 1 / (1 + RMSD_lig_A / 10),     # 0-1, mayor RMSD = menor score
    Score_HB    = ifelse(is.na(HBonds), 0, HBonds / max(HBonds, na.rm=TRUE)),
    Score_total = round((Score_RMSD * 0.6 + Score_HB * 0.4), 3)
  ) %>%
  arrange(desc(Score_total))

candidatos %>%
  select(Label, Ligand, RMSD_lig_A, HBonds, Score_RMSD, Score_HB, Score_total) %>%
  kable(
    col.names = c("Sistema","Ligando","RMSD lig (Å)","H-bonds","Score RMSD","Score HB","Score Total"),
    caption   = "Ranking de candidatos (Score = 0.6×RMSD_inv + 0.4×HBonds_norm)"
  ) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"), full_width = FALSE) %>%
  row_spec(1, bold = TRUE, color = "white", background = "#2A9D8F")
Ranking de candidatos (Score = 0.6×RMSD_inv + 0.4×HBonds_norm)
Sistema Ligando RMSD lig (Å) H-bonds Score RMSD Score HB Score Total
Site III — Rimegepant Rimegepant 1.75 0.864 0.8510638 0.6366986 0.765
Site V — Vibegron Vibegron 6.89 1.357 0.5920663 1.0000000 0.755
Site II — Chlorhexidine Chlorhexidine 3.20 0.000 0.7575758 0.0000000 0.455
Site IV — Atovaquone Atovaquone 7.26 0.261 0.5793743 0.1923360 0.425
Site O — Rimegepant Rimegepant 69.86 0.018 0.1252191 0.0132646 0.080 
# Gráfico de ranking
ggplot(candidatos, aes(x = reorder(Label, Score_total), y = Score_total, fill = Label)) +
  geom_col(width = 0.6, show.legend = FALSE) +
  geom_text(aes(label = Score_total), hjust = -0.2, size = 4) +
  scale_fill_manual(values = site_colors) +
  coord_flip(ylim = c(0, max(candidatos$Score_total) * 1.2)) +
  labs(title    = "Ranking de candidatos farmacológicos",
       subtitle = "Score compuesto: 60% estabilidad (RMSD inverso) + 40% H-bonds",
       x = NULL, y = "Score (0–1)") +
  theme_bw(base_size = 12)

16. Métodos

Preparación de sistemas: Las estructuras de los complejos proteína–ligando se generaron mediante docking molecular con AutoDock Vina sobre la proteína de fusión RSV-F (PDB: 5EA3, residuos 1–513). Los 5 sitios de unión se seleccionaron de la biblioteca DrugCentral (2,331 → 1,020 compuestos FDA-aprobados tras filtrado). Las topologías de los ligandos se generaron con ACPYPE v2023.10.27 (campo de fuerza GAFF2, cargas AM1-BCC).

Dinámica molecular: Proteína parametrizada con AMBER99sb-ildn; solvente TIP3P; caja cúbica con borde de 1.2 nm. Minimización de energía (steepest descent, 50,000 pasos), equilibración NVT (100 ps, 310 K, v-rescale), producción 10 ns (dt = 2 fs). Todas las simulaciones en GROMACS 2025.3 con GPU NVIDIA RTX 5080 (-ntmpi 1 -ntomp 16 -nb gpu -pme gpu -bonded gpu -update gpu).

Análisis: RMSD/RMSF/Rg calculados con gmx rms, gmx rmsf, gmx gyrate. Energía potencial con gmx energy (índice 10). SASA con gmx sasa. Puentes de hidrógeno proteína–ligando con MDAnalysis v2.x (HydrogenBondAnalysis, d(D-A) ≤ 3.5 Å, ángulo D-H-A ≥ 150°). Visualización en R 4.x con ggplot2 + patchwork.

sessionInfo()
## R version 4.3.3 (2024-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 24.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: America/Mexico_City
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] knitr_1.52       rmarkdown_2.32   stringr_1.6.0    reshape2_1.4.5  
##  [5] tibble_3.3.1     kableExtra_1.4.1 patchwork_1.3.2  gridExtra_2.3.1 
##  [9] scales_1.4.0     readr_2.2.0      tidyr_1.3.2      dplyr_1.2.1     
## [13] ggplot2_4.0.3   
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6       jsonlite_2.0.0     compiler_4.3.3     Rcpp_1.1.2        
##  [5] tidyselect_1.2.1   xml2_1.6.0         jquerylib_0.1.4    textshaping_1.0.5 
##  [9] systemfonts_1.3.2  yaml_2.3.12        fastmap_1.2.0      plyr_1.8.9        
## [13] R6_2.6.1           labeling_0.4.3     generics_0.1.4     svglite_2.2.2     
## [17] bslib_0.12.0       pillar_1.11.1      RColorBrewer_1.1-3 tzdb_0.5.0        
## [21] rlang_1.3.0        stringi_1.8.9      cachem_1.1.0       xfun_0.62         
## [25] sass_0.4.10        S7_0.2.2           otel_0.2.0         viridisLite_0.4.3 
## [29] cli_3.6.6          withr_3.0.3        magrittr_2.0.5     digest_0.6.39     
## [33] grid_4.3.3         rstudioapi_0.19.0  hms_1.1.4          lifecycle_1.0.5   
## [37] vctrs_0.7.3        evaluate_1.0.5     glue_1.8.1         farver_2.1.2      
## [41] purrr_1.2.2        tools_4.3.3        pkgconfig_2.0.3    htmltools_0.5.9