Resolution and dose dependence of radiation damage in biomolecular systems

Data from several radiation-damage studies using protein crystals are reanalyzed to determine the local dose- and resolution-dependent decay of diffracted intensity at T = 100 K. The results are inconsistent with long-accepted models and with a proposed linear scaling between maximum tolerable dose and resolution, but are reproduced by a simple, physics-based model.


S1. Fitting of experimental intensity vs resolution and dose data when the crystal orientation was unknown
When the crystal orientation was not given, we assumed that one of the three crystal axes was aligned with the beam direction, and performed three separate fits, assuming alignment of each axis. If the crystal dimensions along the axes were substantially (e.g., an order of magnitude) different, we used the orientation that gave the best fit; the other orientations gave substantially worse fits, regardless of the choice of . When the crystal dimensions along the three axes were comparable (e.g., within a factor of 2), then the best-fit  would change with assumed orientation by at most  0.2, varying with the ratio of the beam FWHM to the crystal size.

S2.1. 2D simulations
For 2D crystal simulations, we used as a unit cell a grey-scale image with a resolution of 1024  512 pixels and with pixel values from 0 to 512 (and a maximum value in the image of ~300). The pixel values correspond to the number of electrons per unit pixel area. Effects of initial static crystal disorder were modeled by randomly rotating the image in each unit cell by 0.5-1.
For each x-ray hit, a small nn (typically 55) pixel interaction region centered at ( , ) ij xy on the crystal is randomly selected and copied into a matrix. To account for the edges of this region, the matrix is extended to a 3n3n matrix by reflecting the original matrix about each edge; e.g., in 1D, an original array "a b c d" becomes "d c b a | a b c d | d c b a". This is the default edge handling method of the scipy library of Python. A Gaussian spatial filter is applied to each of the pixels within the nn interaction region of the extended matrix. The 2D Gaussian kernel applied to a pixel at ( ij x , y ) is given by with  between 0.5 and 1, the small  values causing the Gaussian to decay to near zero in only a few pixels. The intensity of a pixel within the interaction region before and after the kernel is applied to that pixel are r x i ,y j ( ) and where the sum is over all pixels in the extended matrix. The pixel value in the original image matrix is updated after the calculation, while the extended matrix values are not changed. Fig. S7 shows the effect of the Gaussian blur.
Simulations were performed using mm arrays of unit cells with m ranging between 4 and 128, and using unit cells with grid sizes ranging from 6432 to 1024512 pixels. Results with smaller unit cells became largely independent of m for m8. Consequently, m=32 was used in most simulations, with m=16 used for the largest unit cell. Simulations were run until the random hits caused the diffraction peaks in the highest resolution shell to fall below the background level, which corresponded to roughly 5-10 hits/pixel. Each simulation was divided into 30-50 segments, and after each segment the FFT of the crystal was calculated. The intensity of each diffraction peak is calculated by subtracting the squared local background pixel count from the squared peak count. The diffraction peaks were binned into resolution shells defined by upper and lower resolutions. The integrated intensity of the diffraction peaks, normalized by the undamaged intensity, in each resolution shell is plotted against the average number of hits per pixel. Fourier Transform of the initial, undamaged crystal has strong diffraction peaks extending out to the maximum q / resolution of the initial image. After some large number of hits, the real space electron density is blurred throughout the sample, and the Fourier transform decays much more rapidly with q.
Movies S1 and S2 show the evolution of the crystal's electron density and the squared amplitude of its FFT, proportional to the diffracted intensity, with dose.

S2.2. 3D simulations
Computationally much more intensive simulations of radiation damage to 3D crystals were performed to assess whether the whether the shapes of the dose curves   , I q D and the exponent in Eq. 4 depended on dimension. PDB entry 3E4H, for tetragonal crystals of the 29 residue plant protein Cyclotide varv F at 1.8 Å resolution was chosen for the simulation. The 3D electron density map of the asymmetric unit was calculated from the PDB file (with hetero atoms removed) using phenix.fmodel, which created an output in .mtz format. This map was discretized using double-precision floating-point format and between 16 and 128 points in each dimension (limited by the available memory for the simulation) by reading from the .mtz file using the Phenix tool phenix.map_value_at_point, and is shown in Fig. S4. The simulations were performed using the m2.2xlarge cluster instance of Cornell Advanced Computing's Red Cloud, which has 28 cores and 192 GB of RAM. This allowed simulation of crystals with 16168 unit cells when using 128 3 voxels per cell; 161616 unit cells were used for all other simulations. The largest simulations took about 7 hours.

Figure S1
The radiation-damage model of Blake & Phillips as modified by Hendrickson (Blake et al., 1962;Hendrickson, 1976). An initially undamaged crystal region may become either disordered, corresponding to an increased average B factor, or amorphous, corresponding to a complete loss of Bragg scattering. Undamaged crystal becomes disordered with a "rate constant" (fraction per unit dose) 1 k and this disordered crystal becomes fully amorphous with a rate constant 2 k . Undamaged crystal can also directly proceed to the amorphous state at rate 3 k . Previous fits to experimental data typically found 3 0 k  (Hendrickson, 1976;Sliz et al., 2003;Warkentin & Thorne, 2010).  al. (Liebschner et al., 2015) in their Figure 1, reproduced here for convenience.

Figure S4
Experimental data (solid circles) for integrated intensity in resolution shells versus dose for thaumatin crystals at 100 K, measured by Liebschner et al. (Fig. 4 in the original manuscript) (Liebschner et al., 2015). Figure 2 shows the same data normalized in each resolution shell by the zerodose intensity. The solid lines indicate results from simulations assuming (a) a perfect top-hat incident x-ray beam profile and an exponent 1   in Eq. 4; (b) the measured beam profile (Fig. 1    This was used as the unit cell for simulations. The actual crystal of 3E4H has a cubic unit cell of size 82 Å and contains 48 copies of the ASU, and was too large for our simulations.

Figure S7
Gaussian blur with =0.6 applied to a 55 pixel region, with edge handling performed as described in Section S1. Pixel counts correspond to electrons/pixel or electron density. The sum of pixel counts -corresponding to the total number of electrons -is conserved in the blurring. The Gaussian function vanishes beyond 3 from the central pixel.

Figure S8
Results of 2D simulations of radiation damage as described in Section S1 using (

Movie Captions
Movie S1 Video showing the evolution of the electron density of the 2D crystal in Fig. 5 as the number of "hits" increases, generated using our model of radiation damage as a series of Gaussian blurs applied at random locations. Audio as in Video S2.

Movie S2
The evolution of the square of the amplitude of the FFT of the electron density shown in Video S1, corresponding to the evolution of the crystal's diffraction pattern. The audio was generated by mapping the magnitude of the q vector of each FFT peak onto a frequency between 50 Hz (for q=0) and 20 kHz, and generating tones for each with an amplitude proportional to the square of the FFT peak amplitude. The large q peaks / high frequency tones disappear rapidly, while the small q peaks / low frequency tones fade out very slowly.