Towards RIP using free-electron laser SFX data

Here, it is shown that simulated native serial femtosecond crystallography (SFX) cathepsin B data can be phased by rapid ionization of sulfur atoms. Utilizing standard software adopted for radiation-damage-induced phasing (RIP), the effects on both substructure determination and phasing of the number of collected patterns and fluences are explored for experimental conditions already available at current free-electron laser facilities.

1. Introduction

The properties of X-ray free-electron lasers (FELs) and associated sample environments have catalyzed the development of new data collection methods, such as serial femto­second crystallography (SFX). In SFX, a complete dataset can be obtained by exposing thousands of randomly oriented crystals to the X-ray beam (Chapman et al., 2011; Boutet et al., 2012; Redecke et al., 2013; Johansson et al., 2012). One very attractive feature of this method is that crystals of very small sizes are sufficient to obtain high-quality data, allowing data to be collected from smaller crystals. These experiments also offer the exciting possibility of determining structures that are free of the normal types of radiation damage, since the diffracting pulse ends before most of the processes induced by X-ray absorption can occur. To date, SFX crystal structures show no signs of the kinds of radiation damage previously observed at synchrotron sources (Boutet et al., 2012), even though the dose at which the single crystal is exposed can exceed 100 times the dose limit at cryogenic temperature (Owen et al., 2006). At such doses there is enough X-ray fluence to cause ionization of every atom in the sample (Chapman et al., 2011) and the exposure time is about a million times faster than neutralizing electron recapture processes. Despite the lack of secondary radiation damage in SFX data, ionization of atoms in the sample is unavoidable since it occurs within the femtosecond timescale of the pulse (Hau-Riege et al., 2004; Quiney & Nugent, 2011). This change in the electronic configuration of atoms leads to the modification of the atomic scattering factors during diffraction, especially at photon energies near absorption edges (Son et al., 2011, 2011a). It has been theorized that, for X-ray fluences already achievable with existing FELs, this ionization could significantly alter the scattering factors of the heavy atoms, hindering the direct application of anomalous phasing methods (Son et al., 2011, 2013). This is analogous to the situation in synchrotron macromolecular crystallography, where radiation damage, if unaccounted for, can dramatically hinder anomalous phasing (Ravelli et al., 2005). In the SFX case, ionization introduces the exciting possibility of determining phases de novo by varying the fluence and hence the scattering factors of the heavy atoms. This new approach could represent a powerful method of experimental phasing which would not require the modification of the native protein crystals if the sulfur atoms already present in the structure could be used. Since ionization occurs at all sulfur atoms and not just those that are found in disulfide bonds, this kind of electronic modification could be of broader applicability than synchrotron-radiation-induced phasing (RIP), which requires disulfide bonds and/or metals. To avoid confusion with phasing by synchrotron radiation damage at cryotemperatures (Ravelli et al., 2003; Evans et al., 2003), and to highlight the central role of the X-ray intensity, we call this new technique `high-intensity radiation-induced phasing' (HI-RIP).

In this paper it is demonstrated that HI-RIP can in theory be used to determine substructures of radiation damage as well as determine high-quality phases. We show this by simulating two SFX experiments, at high and low X-ray fluence, under experimental conditions that mimic those available at the Linac Coherent Light Source (LCLS) (Emma et al., 2010; Boutet & Williams, 2010). Furthermore, the effect of different experimental parameters such as variable fluences and the number of patterns on the quality of the substructure solution and phasing is explored.

2. Models and methods

To simulate SFX experiments, it is necessary to take into account ionization dynamics and the accompanying changes of atomic scattering factors (Son et al., 2011). In light atoms, after absorbing an X-ray, an inner-shell electron is ionized and the resulting vacancy is relaxed via Auger decay, typically within a few femtoseconds, or fluorescence. A deep inner-shell vacancy in a heavy atom that has more than two subshells undergoes a cascade of Auger decays, ejecting several electrons. At high X-ray intensity, photoionizations may occur sequentially even to the same atom, resulting in highly charged ions. Even though X-ray absorption locally occurs to the single atom, charge states are redistributed to neighboring atoms during or after ionization processes in the femtosecond regime. Electrons ejected from photoionizations or Auger decays induce a cascade of collisional ionization that occurs within about 100 fs. These X-ray ionization dynamics happen stochastically to every single atom in the sample. Once charged states are created, the atomic scattering factors are altered. Assuming that individual atoms are ionized independently, the scattering intensity is evaluated by (see Appendix A for details)

where is the momentum transfer, is the X-ray fluence, and ω is the photon energy. Here, X is one of the atom species in the sample. is the structure factor for X, including neither dispersion correction nor ionization, and is the phase difference between the structure factors for X and Y. In contrast to Son et al. (2011a), here the effects of ionization dynamics and anomalous scattering are considered for all atoms in the sample, which are included in the anomalous coefficients. The atom-specific coefficients of a_X and , and the biatom-specific coefficients of B_XY and C_XY are defined in Appendix A. Note that these coefficients depend on both the photon energy ω and the fluence . Since the dynamical fluctuations of the atomic form factors during the intense X-ray pulse are dramatic (Son et al., 2013), time-averaged form factors are not sufficient to describe the scattering patterns at high X-ray intensity.

The XATOM toolkit (Son & Santra, 2012) provides a way to simulate the electronic damage dynamics and to compute the anomalous coefficients present in equation (1), for various atomic species at a chosen X-ray fluence and energy. In Fig. 1, the fluence dependence of the scattering strengths (coefficient ) for the atomic species of a native protein (i.e. containing only H, C, O, N and S) is plotted, in the case of a 6 keV, 10 fs flat-top pulse. The plot shows that, at low photon fluence, the relative scattering strength is close to unity, which means that the scattering strength is very similar to the `synchrotron' case where the ionization is negligible. No dependence of coefficient on the fluence implies that contribution to the scattering intensity is linearly proportional to the fluence if anomalous coefficients B and C are zero. By contrast, at higher fluence, the scattering strength of S in particular is reduced compared with lighter atoms, owing to its large photoabsorption cross section. Up to about 10¹² photons µm⁻² for example, the high-intensity modification of the scattering factors of the light atoms is still negligible, while the sulfur has lost almost 40% of its effective scattering strength.

Figure 1 Relative effective scattering strength as a function of the X-ray fluence, for different atomic species, at 6 keV.

In principle, it is possible to apply the same strategy as used in the high-intensity multiwavelength anomalous diffraction (MAD) phasing technique (Son et al., 2011). Assuming that all anomalous coefficients are determined theoretically or experimentally, the number of unknowns in equation (1) are 2N [ and for X = ], where N is the number of different atomic species in the sample. To solve these unknowns, a set of 2N linearly independent equations is needed. Since the anomalous coefficients vary not only with wavelength but also with fluences as shown in Fig. 1, one has to measure 2N complete datasets at different fluences. In this sense, it should be called a multi-fluence anomalous diffraction method. However, as discussed in detail later on, it might be difficult to accurately measure the molecular structure factors at several distinct fluences using current SFX technology. Instead, one can simply exploit the bleaching effect of the heavy atoms at two fluence points: low fluence for the un­damaged case, and high fluence for the damaged case. This bleaching effect can be thought of as being analogous to the specific radiation damage observed at synchrotron sources and, as a result, RIP or SIR (single isomorphous replacement) are the natural phasing methods to use. By performing two SFX experiments at different X-ray fluences, it should therefore be possible to determine the heavy-atom positions and then determine phases.

To retrieve the sulfur coordinates from the change of ionization, two X-ray fluences must be chosen: a `low dose' fluence and a high fluence (which may be the highest achievable at the beamline) that modifies the scattering properties of sulfurs while leaving light elements relatively unchanged. Currently, the coherent X-ray imaging (CXI) endstation at LCLS can provide a SFX setup with a beam that can be focused nominally to 0.2 µm of diameter using Kirkpatrick–Baez mirrors, providing up to 0.5 mJ per pulse in the photon energy range 6–10 keV. 6 keV is chosen because the photoabsorption cross section of sulfur is the highest within that energy range, and the maximum fluence per pulse achievable at the focal spot theoretically reaches 10¹³ photons µm⁻². Low fluence pulses can be obtained by attenuating the beam with Si attenuators. Based on the ionization dynamics calculations, here the amount of damage at 10¹¹ photons µm⁻² is judged to be negligible, and this value is set as the lowest intensity.

3. Simulations

The Trypanosoma brucei cathepsin B structure recently solved at LCLS [Protein Data Bank (PDB) code 4hwy ] (Redecke et al., 2013), space group P4₂2₁2, a = b = 125.4 Å, c = 54.56 Å, was used to test the method. This protein consists of 340 residues, and includes 19 S atoms (in 5 methionines and 14 cysteines). The complex-valued structure factors of equation (1) [and hence their phases ] were computed separately for each atomic species X (excluding hydrogen) present in the protein structure using the sfall program (Ten Eyck, 1977; Winn, 2011); Friedel mates were averaged using an ad hoc script. The atom- and biatom-specific coefficients were calculated with the XATOM toolkit and the final scattered intensity for each Bragg reflection was calculated using equation (1). In other words, structure factors without any fluence effects were calculated with sfall, and modified to include high-fluence effects with XATOM. The scattering contribution of the bulk solvent region was generated starting from the PDB using the ano_sfall.com script (Holton et al., 2014) and was summed with the scattering factors of the protein.

The SFX experiment was simulated with partial_sim, which is part of the CrystFEL suite (White et al., 2012, 2013). The program takes a list of fully integrated reflection intensities and generates partial reflection intensities of randomly oriented crystals, adding noise to simulate other measurement errors. The detector geometry was chosen to reproduce the 64 tiles of the Cornell–SLAC Pixel Array Detector (CSPAD) (Hart et al., 2012) installed at the CXI endstation. The sample-to-detector distance was set to 11 cm, giving a resolution limit of about 2.7 Å at the corners of the detector for 6 keV energy. The simulated noise and the beam parameters were selected in order to produce R_split (the error metric commonly used for SFX experiments) (White et al., 2012) and 〈I/σ〉 values as close to a real experiment as possible (see, for example, Redecke et al., 2013). Two sets of 3 × 10⁵ patterns were generated, one with 10¹³ photons µm⁻² (as `high fluence') and another using 10<sup>11</sup> photons µm<sup>−2</sup> (as `low fluence'). The simulated noisy partial intensities were then merged using process_hkl (another program from the CrystFEL suite), using a Monte Carlo integration method (Kirian et al., 2010). To test the convergence of the Monte Carlo procedure, smaller subsets of the complete data set ranging from 1000 to 3 × 10⁵ patterns were processed. The RIP workflow (Nanao et al., 2005) was then utilized to determine heavy-atom substructures and determine phases.

4. Phasing

The first step in experimental phasing by HI-RIP is the determination of the `damage substructure'. RIP substructures are generally qualitatively different from those of anomalous dispersion and isomorphous replacement, and require specialized methods (Nanao et al., 2005; de Sanctis & Nanao, 2012). Briefly, because RIP substructures consist of many weak sites, iterative improvement of substructures and the down-weighting of the damaged (high fluence in HI-RIP) dataset are necessary. The down-weighting of the damaged dataset is carried out in a simple manner in which an initial scale is determined using conventional means, and then down-weighted by the scale factor k. In this approach, a k of 1 is the same as not applying down-weighting. In some cases, both techniques are required for the successful determination of the damage substructure, whereas in others only one of these methods is sufficient. Furthermore, it has been observed that a peak in the correlation coefficient (CC) of SHELXD (Sheldrick, 2010) substructures as a function of scale factor (k) is an excellent predictor of RIP signal, which complements the conventional indicators of signal strength such as R<sub>isomorphous</sub> (Nanao et al., 2005). In comparing the high- and low-fluence datasets, such a relationship is observed (Fig. 2). Because the damage substructure and a reference phases set from a refined model are known, it could also be demonstrated that the F<sub>lowfluence</sub> − F<sub>highfluence</sub> map calculated with model phases by ANODE (Thorn & Sheldrick, 2011) showed peaks, and that these peak heights (called `RIP peak heights') were quite high (up to 25.6σ) over sulfur positions in methionine and cysteines, as expected from Fig. 1. The top ten strongest peaks were over the sulfur positions of residues C107, C219, M131, C154, C122, C136, M138, C119, C158 and C192 with peak heights of 25.6, 25.3, 24.9, 24.1, 28.9, 23.7, 23.2, 22.9, 22.4 and 22.2σ, respectively. Model phases are obviously not available when trying to determine a new structure. Therefore, more practically useful measures of whether there is adequate signal are the substructures quality metrics that are produced by SHELXD. The most frequently used metric in substructure solution is the ratio of CC(all) to CC(weak). In a plot of CC(all) against CC(weak), a contrast between trials with high values of CC(all)/CC(weak) typically indicates that substructures are at least partially correct. Plots of CC(all)/CC(weak) for the simulated data do indeed reveal a contrast in solutions (as shown in Fig. 2), with the best occurring at a k of slightly less than 1.0, which is consistent with that observed in synchrotron RIP (Nanao et al., 2005). Because in this test the `correct' substructure is known from a peak search of the RIP peak map from ANODE, the resultant substructures could also be compared with this reference structure with iotbx.emma (Adams et al., 2010), which revealed that substructures were largely correct at a variety of values of k, with a maximum correctness at k = 0.96 (Fig. 3). Taken together, these results indicate that there is adequate signal in a HI-RIP experiment to determine correct radiation damage substructures, in the same manner as in RIP.

Figure 2 SHELXD correlation coefficient ratios (CCall/CCweak) as a function of k for high fluence = 1 × 10^13 photons µm^−2, low fluence = 1 × 10^11 photons µm^−2, 100000 patterns.

Figure 3 Correctness of SHELXD substructures. The substructure with the highest CCall/CCweak was compared with a reference substructure from ANODE using phenix.emma. Purple sites indicate a correctly identified atom. Black sites indicate that a site was not identified. X and Y are the fractional unit-cell coordinates of the sites, as a fraction of the X and Y axis. The diameter of each site represents the RIP peak height, and is normalized to the most intense peak height in the difference map.

Next, it was determined whether, together with the experimentally determined substructures, the differences between high- and low-fluence datasets were enough to determine phases. The best substructures from each k were input into SHELXE (Sheldrick, 2010) for phase calculation, phase improvement and model building. Initial figure-of-merit weighted phase errors (after one round of phase improvement by SHELXE) were very good at 38°, and this improved and converged to 28° after an additional round of substructure and phase improvement. Unlike in synchrotron X-ray and UV RIP, where the main benefit from substructure improvement is the identification of sites with `negative' electron density in RIP difference maps caused by, for example, side chain rearrangements, these rearrangements are not anticipated on the time scale of an SFX experiment. Therefore, in this case, substructure improvement serves the purpose of identifying weaker sites in the F_{low fluence} − F_{high fluence} map that were not identifiable initially by SHELXD.

Having established that, in the initial testing conditions used here, standard RIP analysis could be used to both determine radiation damage substructures and to produce interpretable electron density maps, the effects of various changes to the experimental conditions were studied. Since the number of patterns that can be collected during an X-ray FEL experiment is often limited due to practical reasons such as limited beam time and/or sample, and the fact that 10¹³ photons µm⁻² is a current upper limit for the photon density in a single FEL pulse, the effect of both parameters on phasing were simulated. It was found that substructure solution and correct phases could be achieved with high fluences of 5 × 10¹² photons µm⁻² even down to 20000 patterns [see, for example, the purple curve in the weighted mean phase error (wMPE) of Fig. 4], a value that is comparable with the average number of patterns required for solving SFX structures with standard methods (Boutet et al., 2012; Barends et al., 2014). With a very large number of patterns (n = 3 × 10⁵), slightly lower fluences of 1 × 10¹² photons µm⁻² could be used as well. It is important to note that the critical parameter is not the ratio of fluences but rather the difference between fluences. This is because the effective scattering strength of sulfur does not decrease linearly with increasing fluence, but its relative change is the highest between 5 × 10¹² and 10¹³ photons µm⁻² (see Fig. 1). In general, a larger number of patterns improves the quality of the HI-RIP solution because of improved averaging of errors in the Monte Carlo integration of intensities.

Figure 4 RIP peak height, best weighted mean phase error (wMPE), and the number of substructures having a correctness greater than 80%, as a function of the number of patterns used. The corresponding low fluence is written in the grey boxes, while the high fluence is shown with different colours (see legend in the bottom-left corner).

5. Discussion and conclusions

The success of the HI-RIP technique relies on the relative strength of the scattering factors of the heavy atoms at the two fluences used, and on the accuracy of the structure factor measurement. The former is mainly related to the specifications of the experimental facility, such as the available photon flux and bandwidth, its shot-to-shot variation and the spatial photon distribution inside the beam profile. The latter is influenced by errors introduced by the SFX technique, such as reflection partiality, inhomogeneous crystal size or quality and non-isomorphism. The Monte Carlo integration of intensities can help to average out these error sources, but a large number of observations are required. Noise was added to the simulations reported here, which yielded data with merging statistics similar to those found with experimental data, but the simulations make several assumptions that might be difficult to achieve using current FEL technology. First, it is possible that the shot-to-shot variation in intensities might contaminate the high-fluence diffraction data with data recorded at lower fluences, reducing the ionization contrast. This could be avoided by equipping the beamline with a beam diagnostic monitor after the interaction region, which could be used to sort the diffraction snapshots as a function of the photon fluence. Second, the actual beam profile is not a top-hat function but has wings of lower intensity surrounding the focus. Crystals which pass through these regions would be exposed to lower intensities than those passing through the focal spot. Currently there do not exist published beam profiles for either of the optical layouts at CXI. We have, however, used our workflow to make crude estimates of the effect of a more Gaussian beam profile, by considering high-fluence datasets created by averaging intensities from snapshots at three different photon flux densities: 10¹³, 5 × 10¹² and 10¹² photons µm⁻². Since the low-fluence data are only affected by a negligible ionization effect, a single set at 10¹¹ photons µm⁻² was used. It was seen that, using 50000 patterns, phasing is still possible as long as less than 40% of the reflection intensity measurements are generated from the tail of the beam. In particular, a dataset composed of 20000 patterns at highest fluence, 20000 at medium and 10000 at low fluence showed a RIP peak height of about 21 and a best wMPE of about 20° (roughly twice as large as in the pure beam case). Doubling the number of low-fluence patterns at the expense of higher fluence data reduces the RIP contrast much below 16σ, so that no phasing is possible. This indicates that the HI-RIP method can tolerate a more Gaussian beam profile. The simulation, however, does not fully model the case of a Gaussian beam of the same (or smaller) size as the crystalline sample, where the single diffraction pattern will contain both high-fluence and low-fluence scattering. This condition must be modeled with a new set of equations, by introducing different electronic responses to various X-ray fluences over the sample. Experimentally, these problems could potentially be avoided by using an X-ray beam size much larger than the crystal size. To treat the ionization dynamics in the present work, an independent atomic model, ignoring charge rearrangement with neighboring atoms, has been used. The molecular environment will affect all anomalous coefficients calculated within the independent atomic model. In our model, resonant absorption processes, shakeup or shakeoff processes (Persson et al., 2001), and collisional ionization (Caleman et al., 2009), which induce further ionization are also neglected. These assumptions, and no doubt others, might cause a discrepancy between simulated data and experimental data and will be incorporated into future simulations.

To conclude, SFX data at high and low photon fluence have been simulated and it has been shown that standard RIP software can be used to both determine substructures and produce interpretable electron density maps. Furthermore, it has been demonstrated that phasing can be successful using wide ranges of fluences and numbers of patterns. Future work should include an experimental demonstration of the technique, taking into account the above points and using the simulations presented here as a guide for planning.

APPENDIX A

Scattering intensity including ionization dynamics for all atomic species

Including all types of ionization dynamics, the scattering intensity (per unit solid angle) can be evaluated by (Son et al., 2011, 2013)

where g(t) is the normalized pulse envelope, is a coefficient given by the polarization of the X-ray pulse, and P_I is the population of the global electronic configuration I. This global configuration index is given by I = , where X is one of the atom species in the sample and I_j ^X indicates the electronic configuration of the jth atom of X, which is located at the position . Assuming that individual atoms are ionized independently, the population of configuration I can be decomposed into individual populations, = . Using this assumption, one can expand (2) and group terms for the same species, which yields equation (1). In (1), the coefficients of a_X, , B_XY and C_XY are a function of , defined by equations (3),

where I_X is the electronic configuration of the atomic species X and indicates the normal atomic form factor of the ground configuration of X, which includes neither dispersion correction nor ionization. = × is the atomic form factor of X at given time t, describing synchronized changes in the atomic form factor over all elements of the same kind. Among these coefficients, is corresponding to the relative effective scattering strength of X, as plotted in Fig. 1.