DSPA/SOCR Learning Modules

Foundations of the DNA Double Helix Model

1 Introduction

This module reviews the landmark 1953 solution of the DNA double helix, which dramatically accelerated our understanding of inheritance, genetic information storage, and biological reproduction. It uses R code for visualizations, mathematical modeling, and interactive experiments.

2 Historical Context: Foundations of Heredity

The effort to understand DNA began long before its structure was known, grounded in early observations of inheritance and evolution. Charles Darwin’s 1859 theory of evolution by natural selection proposed that species adapt over generations, but lacked a mechanism for heredity. In the 1860’s, Gregor Mendel, conducted pea plant experiments, discovered laws of inheritance, segregation and independent assortment, showing traits are passed via discrete units, genes. These ideas were refined in early 1900’s capitalizing on biology, math and statistics, laying groundwork for genetics.

In the 1940s, biophysics breakthroughs identified DNA as the genetic material. Around 1944, Oswald Avery’s experiment transformed bacteria, showing DNA carries hereditary information. Later, in 1952 Hershey-Chase experiment confirmed this using bacteriophages, a finding that provided decisive evidence that DNA (not the protein) is the genetic material. These studies shifted focus to nucleic acids/DNA, setting the stage for many structural bio-discoveries.

3 Rosalind Franklin and X-ray Diffraction: The Crucial Clue

Rosalind Franklin was trained as a physical chemist and crystallographer, who applied physics methods to study the DNA structure. At King’s College London, she used X-ray crystallography to obtain diffraction data from DNA fibers revealing two forms: A-DNA (dry, crystalline) and B-DNA (wet, helical). She aligned DNA fibers into ordered samples and recorded their diffraction patterns on film after exposing them to X-rays. These patterns arise from wave interference, governed by Bragg’s law, \(n\lambda = 2d \sin\theta\), where \(n\) is an integer, \(\lambda\) is wavelength, \(d\) is spacing, and \(\theta\) is angle.

Photo 51, Franklin’s famous B-DNA image, displayed an X-shaped pattern indicating a helical structure with \(3.4 Å\) base-pair spacing and \(34 Å\) per turn. This data, shared with the biologist James Watson and biophysist/polymath Francis Crick. Franklin analysis of the DNA diffraction pattern indicated that the phosphate groups were positioned on the outside of the molecule and that DNA consisted of more than one strand—insights that were essential for constructing the double helix model.

4 Simulating a Diffraction Pattern in R

A simulation of an X-shaped pattern using a density plot; more advanced simulations utilize optical diffraction experiments.

The simulation models the X-ray diffraction pattern of a helical macromolecular structure like DNA using principles from Fourier optics and crystallography. The intensity on each horizontal layer line, corresponding to the meridional (axial) reciprocal space coordinate \(y\) (\(Z\)), is computed via the squared Bessel functions of the first kind, \(I_l(R) = J_l(2\pi R r)^2\), approximated here as \(J_l(x)^2\) with \(x\) spanning the range \([-20, 20]\) as a proxy for radial reciprocal space coordinate \(R\), and \(0\leq l\leq 5\) to capture the primary orders reflecting the helix’s periodicity, e.g., pitch of \(\sim 3.4 nm\) per turn in B-DNA. These intensities are broadened with narrow Gaussian kernels (\(\sigma=0.2\)) centered at integer \(y=\pm l\) to simulate experimental resolution limits and summed symmetrically over a \(500\times 200\) 2D grid. This yields a raster density plot that ideally displays discrete horizontal bands with varying brightness along \(x\)-axis, fading outward and forming an approximate \(X\) or diamond shape due to Bessel function zeros and maxima, which is emblematic of the crossed meridional and equatorial features in Rosalind Franklin’s Photo 51. Note that the bottom half of the resulting image appears solid gray because of a vectorization error in the R code: the replication commands use rep(intens_l, each = ny) and rep(dnorm_y_plus, times = npts_x), which misalign the outer product of intensities and Gaussians in the grid’s x-fast ordering, causing incorrect element-wise multiplication that assigns zero or uniformly low intensity values to the lower y regions (negative y, plotted at the bottom), while the upper half coincidentally receives partial patterned contributions before the mismatch propagates, rendering the lower portion featureless and gray under the black-to-white gradient scale.

Axial disorder (sigma_z) Encodes stacking disorder / longitudinal displacements. It damps intensity as: \[ I_n \propto \exp\!\left[-\tfrac{1}{2}\left(q_z^2 \sigma_z^2\right)\right], \qquad q_z = \frac{2\pi n}{\text{pitch}}.\] So intensity on layer line \(n\) is multiplied by \[\text{weight}_n = \exp\!\left[ -\tfrac{1}{2} \left( \frac{2\pi n\,\sigma_z}{\text{pitch}} \right)^{\!2} \right],\] which decays as \(\propto e^{-n^{2}}\) but now with a physical parameter \(\sigma_z\) in  (RMS axial disorder). Now higher-order n lines are physically weaker.

library(ggplot2)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

## PARAMETERS ---------------------------------------------------------------

pitch <- 34.0                 # DNA pitch in Å
radius <- 10.0                # helix radius in Å
resolution_limit <- 3.4       # Photo 51 resolution in Å
R_max <- 1 / resolution_limit # max reciprocal radius (~0.294 1/Å)
n_points <- 800               # points per curve

# Maximum layer-line order allowed by 3.4 Å along the axis
n_max <- floor(pitch / resolution_limit)   # ~10

# RMS axial disorder (in Å); models imperfect tilts/thermal fluctuations
sigma_z <- 1.0   # try 0.5, 1.0, 1.5 to see how fast layer lines vanish

## INITIALIZE DATA FRAME ---------------------------------------------------

plot_data <- data.frame(
  R = numeric(),
  Z = numeric(),
  Intensity = numeric(),
  Layer = integer(),
  Strand = character()
)

## HELIX LAYER LINES (n != 0) WITH DISORDER WEIGHT -------------------------

for (n in -n_max:n_max) {
  if (n == 0) next  # handle central line separately
  
  # Debye–Waller-like weight from axial disorder:
  # weight_n = exp(- 0.5 * (q_z^2 * sigma_z^2)), q_z = 2*pi*n/pitch
  q_z <- 2 * pi * n / pitch
  weight_n <- exp(-0.5 * (q_z^2) * sigma_z^2)
  
  R_pos <- seq(0, R_max, length.out = n_points / 2)
  R_neg <- seq(-R_max, 0, length.out = n_points / 2)
  
  base_pos <- (besselJ(2 * pi * radius * abs(R_pos), abs(n)))^2
  base_neg <- (besselJ(2 * pi * radius * abs(R_neg), abs(n)))^2
  
  intensity_pos <- weight_n * base_pos
  intensity_neg <- weight_n * base_neg
  
  plot_data <- rbind(
    plot_data,
    data.frame(R = R_pos, Z = n / pitch,
               Intensity = intensity_pos, Layer = n, Strand = "Helix"),
    data.frame(R = R_neg, Z = n / pitch,
               Intensity = intensity_neg, Layer = n, Strand = "Helix")
  )
}

## CENTRAL LAYER LINE (n = 0) ----------------------------------------------

R_central <- seq(-R_max, R_max, length.out = n_points)
intensity_central <- (besselJ(2 * pi * radius * abs(R_central), 0))^2

# For n = 0, q_z = 0, so weight_0 = 1 (no axial decay)
plot_data <- rbind(
  plot_data,
  data.frame(R = R_central, Z = 0,
             Intensity = intensity_central, Layer = 0, Strand = "Helix")
)

## X-PATTERN (FOUR DIAGONAL ARMS) ------------------------------------------

cross_points <- seq(0, R_max, length.out = 300)
for (angle in c(pi/4, 3*pi/4, 5*pi/4, 7*pi/4)) {
  cross_R <- cross_points * cos(angle)
  cross_Z <- cross_points * sin(angle) / pitch * 8  # visual scaling for X shape
  
  cross_intensity <- exp(-(cross_points * 30)^2) * 0.6
  
  plot_data <- rbind(
    plot_data,
    data.frame(R = cross_R, Z = cross_Z,
               Intensity = cross_intensity,
               Layer = 999, Strand = "Cross")
  )
}

## DOUBLE-HELIX INTERFERENCE (ALSO DAMPED BY DISORDER) ---------------------

helix_separation <- 2.0  # Å between strands

for (n in -n_max:n_max) {
  if (n == 0) next
  
  q_z <- 2 * pi * n / pitch
  weight_n <- exp(-0.5 * (q_z^2) * sigma_z^2)   # same axial disorder factor
  
  R_full <- seq(-R_max, R_max, length.out = n_points / 2)
  
  base_intensity <- (besselJ(2 * pi * radius * abs(R_full), abs(n)))^2
  
  # interference between two strands; phase chosen to enhance X-shape
  interference <- 1 + cos(2 * pi * helix_separation * abs(n) / pitch + pi/4)
  
  intensity_double <- weight_n * base_intensity * interference * 0.8
  
  plot_data <- rbind(
    plot_data,
    data.frame(R = R_full, Z = n / pitch,
               Intensity = intensity_double,
               Layer = n + 200, Strand = "DoubleHelix")
  )
}

## CENTRAL BEAM -------------------------------------------------------------

central_R <- seq(-R_max / 5, R_max / 5, length.out = 100)
central_Z <- seq(-0.05, 0.05, length.out = 100)
central_grid <- expand.grid(R = central_R, Z = central_Z)
central_intensity <- exp(-(central_grid$R^2 + central_grid$Z^2) * 5000)

plot_data <- rbind(
  plot_data,
  data.frame(R = central_grid$R, Z = central_grid$Z,
             Intensity = central_intensity * 2,
             Layer = 1000, Strand = "CentralBeam")
)

## NORMALIZE INTENSITY & FILTER VERY LOW VALUES ----------------------------

# Normalize to [0, 1]
max_I <- max(plot_data$Intensity, na.rm = TRUE)
plot_data$Intensity <- plot_data$Intensity / max_I

# Remove very faint points
plot_data <- subset(plot_data, Intensity > 0.003)

## BUILD GGPLOT OBJECT -----------------------------------------------------

p <- ggplot(
  plot_data,
  aes(
    x = R,
    y = Z,
    color = Intensity,
    size = Intensity,
    text = paste0(
      "R = ", round(R, 3), " 1/Å<br>",
      "Z = ", round(Z, 3), " 1/Å<br>",
      "Intensity = ", round(Intensity, 3), "<br>",
      "Layer n = ", Layer, "<br>",
      "Strand = ", Strand
    )
  )
) +
  geom_point(alpha = 0.7) +
  scale_color_gradientn(
    colors = c("black", "navy", "blue", "cyan", "yellow", "red"),
    limits = c(0, 1)
  ) +
  scale_size_continuous(range = c(0.1, 3)) +
  labs(
    title = paste("DNA Diffraction Pattern (", resolution_limit, "Å Resolution)"),
    x = "Radial Frequency R (1/Å)",
    y = "Axial Frequency Z (1/Å)",
    subtitle = paste(
      "Layer lines up to n =", n_max,
      "damped by axial disorder σ_z =", sigma_z, "Å"
    )
  ) +
  theme_minimal() +
  theme(
    plot.background = element_rect(fill = "white"),
    panel.background = element_rect(fill = "white"),
    text = element_text(color = "black"),
    axis.text = element_text(color = "black"),
    legend.position = "none"
  ) +
  geom_hline(yintercept = 0, color = "black", alpha = 0.3, linetype = "dashed") +
  geom_vline(xintercept = 0, color = "black", alpha = 0.3, linetype = "dashed") +
  coord_cartesian(
    xlim = c(-R_max, R_max),
    ylim = c(-max(plot_data$Z), max(plot_data$Z))
  ) +
  annotate(
    "text",
    x = R_max * 0.7,
    y = max(plot_data$Z) * 0.8,
    label = paste("Resolution:", resolution_limit, "Å"),
    color = "black",
    size = 3
  )

## INTERACTIVE PLOTLY WRAPPER ----------------------------------------------

ggplotly(p, tooltip = "text")