Research ArticleGEOPHYSICS

Lateral variation of the Main Himalayan Thrust controls the rupture length of the 2015 Gorkha earthquake in Nepal

See allHide authors and affiliations

Science Advances  26 Jun 2019:
Vol. 5, no. 6, eaav0723
DOI: 10.1126/sciadv.aav0723


The Himalaya orogenic belt produces frequent large earthquakes that affect population centers along a length of over 2500 km. The 2015 Gorkha, Nepal earthquake (Mw 7.8) ruptured the Main Himalayan Thrust (MHT) and allows direct measurements of the behavior of the continental collision zone. We study the MHT using seismic waveforms recorded by local stations that completely cover the aftershock zone. The MHT exhibits clear lateral variation along geologic strike, with the Lesser Himalayan ramp having moderate dip on the MHT beneath the mainshock area and a flatter and deeper MHT beneath the eastern end of the aftershock zone. East of the aftershock zone, seismic wave speed increases at MHT depths, perhaps due to subduction of an Indian basement ridge. A similar magnitude wave speed change occurs at the western end of the aftershock zone. These gross morphological structures of the MHT controlled the rupture length of the Gorkha earthquake.


The collision of the Indian and Eurasian plates formed the Himalaya, the highest mountains in the world, and presents a serious earthquake hazard to millions of people (1). Many important science questions remain about Himalayan formation and seismic hazards, including distinguishing between different possible geometries of the Main Himalayan Thrust (MHT) (2, 3) and better definition of structural causes and locations of rupture segmentation both across-strike and along-strike of the orogenic belt (46). Although previous studies have investigated fault morphology and structural heterogeneity beneath the High Himalaya, the complexity and lack of near-field observations makes interpretations uncertain. The devastating 25 April 2015 Mw (moment magnitude) 7.8 earthquake with its epicenter in Gorkha, Nepal has drawn renewed attention to the severe Himalayan geohazards anticipated from future large earthquakes. Our seismic observations of the aftershock zone provide a unique opportunity to further understand the frontal part of the continental collision zone (Fig. 1).

Fig. 1 Tectonic map with seismicity and seismic stations.

Earthquakes include the Mw 7.8 Gorkha and Mw 7.3 Kodari earthquakes (yellow stars), aftershocks (red circles), earthquakes that occurred before the Gorkha earthquake but after 1990 (gray circles) (, and historic seismicity of Mw >7.5 since 1500 (large solid gray circles). Seismic stations include those we deployed along the China-Nepal border (blue triangles), in the aftershock zone (black triangles), and those of the Nepalese Seismic Network (green triangles). EVN (Everest) is a broadband seismic station deployed at the Khumbu Valley with arrival time and waveform data available since May 2015. The arrow indicates the direction of the Indian plate with respect to the Eurasian plate. The thick gray curve shows the location of Hi-CLIMB (Himalayan-Tibetan Continental Lithosphere During Mountain Building) (3) and HIMNT (Himalaya Nepal-Tibet Seismic Experiment) (8) cross sections (Fig. 4) and stations. The gray rectangle indicates the area shown in Fig. 3. Inset: Seismic stations we used from the larger area. MFT, Main Frontal Thrust; MBT, Main Boundary Thrust; MCT, Main Central Thrust; STD, South Tibet Detachment; and YZS, Yarlung Zangbo Suture.

After the recognition of the MHT based on focal depths and fault plane solutions of moderate earthquakes (2), studies of the MHT proliferated. The MHT has been variously identified on the basis of wave speed changes likely due to water released from underthrust sediments and trapped beneath the fault, lithologic changes (3, 7), an anisotropic zone perhaps due to shearing above the fault (8), intense microseismic activity (9), and high electrical conductivity (10). The source parameters of the Gorkha earthquake given by the common earthquake catalogs such as from the U.S. Geological Survey (USGS) show that the rupture nucleated at ~8-km depth and propagated for ~120 km along-strike and ~80 km up- and down-dip. However, this model appears rather simplistic when compared to the geometry of the MHT interpreted from geophysical and geological studies (7, 912), all of which show a characteristic ramp-flat-ramp thrust-belt structure across-strike.

The mainshock was followed by a large number of aftershocks including the 12 May 2015 Mw 7.3 Kodari aftershock (13) that was initiated near the eastern termination of the mainshock rupture. Although it is widely accepted that the mainshock nucleated and ruptured on the MHT, the aftershock activity has remained controversial. Aftershocks have been relocated above the MHT and treated as an illuminating faulting structure in the Himalayan orogenic prism (1416) or on the MHT, nominally delineating a ramp geometry of the MHT (17), and even below the MHT, interpreted as related to earthquake-induced stress transformation in the subducting Indian lower crust (18, 19).

The 2015 Gorkha earthquake is the largest earthquake to occur on the central segment of the Himalayas in 80 years and the first major event since digital seismic observations have been available. We combine seismic waveforms from several different deployments, including our network of 22 seismic stations deployed in China along the China-Nepal border to the north of the aftershock zone before the mainshock and our NAMASTE (Nepal Array Measuring Aftershock Seismicity Trailing Earthquake) network of 46 stations within and to the south of the aftershock zone 2 months after the mainshock. We also used waveform and arrival time data from the China National Seismic Network, the Nepalese Seismic Network, and other stations in nearby countries (Fig. 1). Using arrival times and waveform modeling (Fig. 2 and fig. S1), we determined source parameters of earthquakes, velocity structures, and discontinuity topography in and around the source area (Fig. 2; see Materials and Methods).

Fig. 2 Data for an example earthquake on 23 August 2015.

Record section of vertical-component seismograms recorded by (A) NAMASTE stations and (B) China-Nepal border stations used to read P- and S-wave arrival times. (C) Waveform inversion for the determination of a fault plane solution. Pz, vertical component of the P wave; Pr, radial component of the P wave; Sz, vertical component of the S wave; Sr, radial component of the S wave; Sh, transverse component of the S wave. (D) Total of 533 earthquakes in the Gorkha earthquake sequence we studied here, which occurred each month since April 2015.

Our results show substantial cross-strike and along-strike variations in the inferred geometry of the MHT and in wave speed variations around the MHT. Our results include many features of previous studies of earthquake relocations (14, 15), fault plane solutions (17), tomographic models (2023), and receiver function images (3, 24), but our new model has better depth and lateral resolution because of our combination of all the datasets mentioned above. We took advantage of the merged dataset by adopting a comprehensive strategy using direct, reflected, and refracted waveforms recorded by local, regional, and teleseismic stations as discussed in Materials and Methods.


We use our 266 well-located earthquakes (table S1), 18 well-determined focal mechanisms (Fig. 3 and table S2), and our three-dimensional tomographic model (data file S1) to interpret the earthquake activity and velocity structure of the Himalayan collision zone, particularly in the depth range of the MHT (Fig. 3). These 266 earthquakes are chosen from the total of 533 events (Mw > 3.5) in the Gorkha earthquake sequence because they have more than 30 recording stations, and relocation uncertainties estimated to be less than ±0.5 km (fig. S2). The epicenters of most of these earthquakes are clustered along the boundary between the physiographic Lesser and Higher Himalaya trending west-northwest (WNW) to east-southeast (ESE), the so-called PT2 (physical transition 2), locally coincident with the Gorkha-Pokhara Anticlinorium (12). Mechanisms for most of the events show reverse faulting, while some of the earthquakes located between the Gorkha earthquake epicenter and the Kathmandu basin exhibit strike-slip faulting (table S2). In the western part of the aftershock zone (longitude <85.5°E), there are two separate seismic belts located south and north of the Kathmandu basin (25). Between these two seismic belts, in the area corresponding to large coseismic slip in the beginning 30 s of the mainshock rupture process (4, 26), aftershocks are rare. In contrast, to the east (longitude >85.5°E), most of the aftershocks are concentrated in a single belt southwest of the Mw 7.3 Kodari earthquake, suggesting along-strike variation in fault structure (fig. S3).

Fig. 3 Map view of our earthquake relocations, fault plane solutions, and P-wave velocity structure for the MHT surface.

The white lines indicate the locations and the direction (25° and 115°, respectively) of the cross sections shown in Fig. 4. Two gray patches (upper left and lower right) show the region of large slips [50% of the maximum slip, as is usually defined as the main slip area for large earthquakes (13)] of the 2015 Mw = 7.8 Gorkha (4) and the 1934 Mw = 8.2 Bihar-Nepal (27) earthquakes. Two wide gray lines are gravity lineaments (30). Black dots are well-relocated earthquakes. The beach balls are the focal mechanisms of earthquakes obtained from waveform modeling. The red square is Kathmandu.

Our new tomographic images on the MHT surface (see Fig. 3 and data file S1) show high Vp values in the aftershock zone, similar to previous tomographic images at shallow depth (23). However, because of many arrival times recorded with dense station coverage in and around the aftershock zone, we are able to image velocity variations in more detail, as confirmed by our resolution tests (fig. S4). The following detailed features are visible: (i) The high-Vp zone (~6.3 km/s) extends ~250 km along-strike, 100 km beyond the aftershock zone to ~87°E where the 1934 Mw = 8.2 earthquake occurred but is bounded to WNW and ESE by lower wave speeds (~5.8 km/s). (ii) There are two cross-strike bands of high-Vp zones at ~84.7°E and ~86.5°E. (iii) The low-Vp region between the two cross-strike high-Vp zones is in the area where the southern aftershock belt is absent.

Cross-strike and along-strike cross sections show the along-strike variations in the Gorkha aftershock zone (Fig. 4). We locate the hypocenter of the Mw 7.8 earthquake ~16.5 km below sea level (we give all depths as below sea level), comparable with the depth of the MHT estimated from previous geophysical (3, 10, 24) and geological (12) studies (Fig. 4A). Both the Mw 7.8 mainshock and the Mw 7.3 largest aftershock have similar north dipping nodal planes dipping ~5°N to 7°N (13) (Figs. 3 and 4), strongly suggestive of rupture of the MHT, although likely an oversimplification as representing an average of MHT dip across the rupture area. Our ~18.5-km rupture depth for the Mw 7.3 aftershock is slightly deeper than the mainshock, in good agreement with previous estimates of earthquake location from the Nepalese Seismic Network (25) and the MHT depth from receiver function analysis (Fig. 4B) (8). Aftershocks become deeper to the east as seen in Profile aa′ for over 100 km (Fig. 4C), above a similar but steeper eastward deepening of the Moho in our tomographic image, consistent with previous geophysical studies (21, 24).

Fig. 4 Cross sections across-strike and along-strike of the collision zone (for profile locations see Fig. 3).

The background in (A) and (B) is receiver function images along Hi-CLIMB and HIMNT (Fig. 1) (3, 8), and in (C) is our P-wave tomography with contour lines (0.5 km/s). Earthquakes are projected no more than 50 km into the section. Black lines show our preferred MHT and Moho. Other symbols show constraints on MHT and Moho from receiver functions (dashed lines): N (3) and H (24); tomographic images (colored solid lines): K (21) and M (20); MT soundings (green broken dotted line): L (10); and geological reconstruction (blue broken dotted line): J (12) and cyan lines above the MHT in (A) (12, 27).


Locating earthquakes in plate convergence regions around the world is problematic due to the uncertainties in the local three-dimensional velocity structure and lack of near-source recordings. Multievent relocation is an improved method to precisely determine relative hypocenter locations. However, the absolute values of the relocated hypocenters can still be affected by the assumed velocity model and starting earthquake locations (15, 18). Earthquake locations could be deeper than we show if the MHT and Moho discontinuities are shallower and inferred P- and S -wave velocities are higher than the actual values. Our careful choice of a suitable velocity model and our use of different types of waveform should yield accurate earthquake source parameters.

As expected, the MHT deepens to the north-northeast (NNE) parallel to the plate convergence direction as the Indian plate subducts beneath Tibet, as does the Moho transition zone, shown as the region between isopleths of Vp = 7.8 and 8.0 km/s (Figs. 4 and 5). Both the MHT and Moho have average dips of ~5° NNE across strike, consistent with the nodal plane of the Gorkha earthquake. It is remarkable, however, that the Moho, as defined by Vp isopleths, also has an average dip of >5° ESE, along-strike, for more than 100 km (Fig. 4C). Although the Moho depth in our tomographic images [and others (20, 21)] is systematically deeper than that obtained by receiver function analysis (Fig. 4) (3, 24), likely due, in large part, to our arbitrary choice of 7.8 to 8.0 km/s as marking the Moho, comparison of these two profiles shows the same ESE dip (Fig. 4C). The deeper edge of the distribution of the aftershocks to the east of 85°E (i.e., excluding the mainshock and a few nearby events) also dips ESE, although only ~2°. This observation is reinforced by a similar ESE dip of the 6 km/s isovelocity line. Because of the nature of our data, we cannot distinguish whether the smooth along-strike dips in our tomography and aftershock catalog (Fig. 4C) may actually conceal lateral steps in the MHT or Moho depth.

Fig. 5 Final velocity model showing map view of the MHT and a cross section passing through the Gorkha mainshock hypocenter.

Front face to the east and rear face to the west of the cube are Vp wave speeds below and above the MHT surface along the line of profile through the Gorkha mainshock at 84.8°E. Colored subhorizontal plane is the wave speed as measured along our interpreted MHT, showing clear along-strike variability. The wide gray line is gravity lineament (30). Dashed line is a 4000-m elevation contour that trends SE cross-cutting the ESE orogenic strike (elevations increase eastward). Image is true scale (depth scaled 1:1 with latitude and longitude).

About 95% of our well-located aftershocks are shallower than the Mw 7.8 and 7.3 events; hence, they are located within the hanging wall above the MHT. Thus, the deeper edge of the main aftershock distribution likely marks the MHT. The great preponderance of earthquakes in the hanging wall, compared to the small number of earthquakes in the footwall, is consistent with neotectonic evolution of the MHT by duplexing that created many faults in the hanging wall capable of reactivation (7). The southern and northern seismic belts in the western aftershock zone have lower edges at ~11.5 and ~16.5 km, respectively, with only rare aftershocks in-between (a to a′ in Fig. 4). Previous studies have mentioned a steep ramp beneath the Lesser Himalaya within the western aftershock zone (3, 10, 12, 24), hereafter called the Lesser Himalayan ramp (LHR), to distinguish it from the Higher Himalayan ramp further north that is likely much larger but deeper so aseismic (27). We show that the LHR is 25 to 30 km wide and dips between 9° and 12° NNE, which is gentler than the ramp suggested by some previous studies (10, 12, 24), although slightly steeper than the fault plane of the Gorkha earthquake (Fig. 4A) (13). The southern seismic belt terminates to the east beneath the Kathmandu allochthon. In contrast, earthquakes in the eastern aftershock zone and around the Mw 7.3 largest aftershock are more concentrated and terminate downward at a rather uniform depth, suggesting a flatter and deeper MHT than to the west (Fig. 4B).

Although abrupt lateral transitions in the geometry of the MHT between the source areas of the 1505 Mw 8.1 and 1934 Mw 8.2 earthquakes have been proposed from structural reconstructions and thermomechanical modeling (11), this geometry has never been directly imaged. Here, we show that the three distinct lateral transitions recognized from earthquake distributions (western: at the Gorkha mainshock and the western limit of rupture; central: at the transition from two seismic belts to one, close to Kathmandu; and eastern: at the eastern limit of rupture close to the largest aftershock) are all spatially related to lateral wave speed variations and hence to along-strike structures.

High Vp values as we find in the aftershock zone are likely associated with high crustal strength that, in turn, can sustain large earthquakes. The Gorkha mainshock hypocenter is located at the western edge of the along-strike high-Vp region, and the 1934 Mw 8.2 Bihar-Nepal earthquake occurred at the eastern end of the same high–wave speed region (28). Thus, the 2015 and 1934 earthquakes all initiated close to steep lateral velocity gradients (Fig. 3A) that likely mark structural boundaries between differing geologic blocks.

The low-Vp region that coincides with the termination of the southern aftershock belt is too deep to correspond to the Kathmandu basin so likely represents either or both a laterally restricted thick Precambrian basin beneath the MHT (29) or fluids trapped beneath the MHT (10). Our Vp/Vs tomography shows that the low-Vp region is also a high-Vp/Vs zone (fig. S5A) and thus more likely to correspond to fluids trapped in thin cracks in crystalline basement. This along-strike wave speed gradient and termination of the southern aftershock belt are likely linked to the cross-strike structure of the MCT that steps >50 km toward the foreland (Fig. 3).

The cross-strike high-Vp anomaly to the east of the aftershock zone has its western edge aligned with a cross-strike basement structure inferred from gravity data (30), the “Patna fault Tangra Yum Co rift” lineament (Figs. 3 and 5). Our high wave speeds may represent the Munger-Saharsa basement ridge directly beneath the MHT, although our observed high Vp does not extend as far east as the “Kishangang fault–Xainza Pum Qu rift” lineament thought to be the eastern boundary of the Munger-Saharsa ridge (29, 30). Reactivation of the Patna faults during Indian plate subduction (22) may provide the lateral ramp marked by our increased wave speeds that limited the eastern propagation of the Gorkha earthquake (12). Many aftershocks in the orogenic prism exhibit steeper dip angles than the MHT, indicating reactivation of duplex structures above the MHT in this context (27). In addition, near-vertical active faults exist in the Kathmandu basin, which may correspond to strike-slip earthquakes at MHT depths. Future seismic imaging should address the true nature of these lateral transitions at the MHT and even Moho depths so that we can better predict the occurrence of other lateral changes that may have the potential to control earthquake rupture patterns. The 2015 Gorkha earthquake has only ruptured a small portion at the eastern end of the nearly 800-km-long seismic gap, which has had no large earthquakes for over 500 years (1). Following our observations of the potential for structural barriers to control major earthquake ruptures, the impetus now must be to image the entire 2500-km-long Himalayan front to determine the morphology of the MHT and the likely controls on the maximum magnitude of rupture that can be accommodated in different parts of this convergent zone.



Our estimates are based on five datasets: (i) waveform data from 22 broadband seismic stations along the China-Nepal border, deployed by the Institute of Tibetan Plateau Research, Chinese Academy of Sciences before the Gorkha earthquake (blue triangles in Fig. 1) (15); (ii) waveform data from a dense seismic array of 45 seismic stations, deployed by the NAMASTE project team 2 months after the mainshock (black triangles in Figure 1) (31); (iii) waveform and arrival time data from the Data Management Centre of China National Seismic Network (SEISDMC; doi:10.11998/SeisDmc/SN) (32) and the China Earthquake Network Center (CENC) (blue triangles in the inset in Fig. 1); (iv) bulletins from the Nepalese Seismic Network for earthquakes that occurred since 1990 when most of the network was constructed (green triangles in Fig. 1) (9, 14); (v) bulletins from the USGS Data Management Centre (DMC) and International Seismological Centre for earthquakes of Mw > 3.5 that occurred since 1990 (black triangles in the inset in Fig. 1). We used arrival time data from 1247 earthquakes, comprising the Gorkha mainshock, 532 aftershocks of Mw > 4.0, and 714 earthquakes that occurred in the whole study area before the Gorkha earthquake since 1990 to study the source parameters and seismic tomography (red and gray circles in Fig. 1). These data include approximately 15,000 P- and 6000 S-wave arrival times within epicentral distances of 1000 km (fig. S1).

Earthquake relocation

We calculated the earthquake hypocenters in two steps. First, we determined hypocenters for the Gorkha earthquake and 532 aftershocks based on the HYPOSAT methodology (33). The combination of local and regional seismic data is optimal for precise hypocenter determination. Using the hypocenters determined by HYPOSAT, we calculated hypocenters of all earthquakes using a multiscale double-difference earthquake relocation method (multi-DD) (34), which is modified from the hypoDD programs (35), to include arrival times of other phases recorded by regional seismic networks. The original hypoDD program assumes a flat Earth model and is appropriate for local-scale calculations. It minimizes residuals between observed and theoretical travel time differences, drkij for pairs of earthquakes i and j observed at the same station k, thusdrkij=(dtkij)obs(dtkij)cal=(tkitkj)obs(tkitkj)cal(1)

For regional scale, the sphericity of Earth should be taken into account (34). For eight moderate aftershocks (Mw 5.5 to 6.7 in table S2), we constrain their focal depths by modeling the teleseismic waveforms of the direct P and the surface reflection pP and sP waves (15, 36, 37). To help constrain the absolute focal depths of other earthquakes, we hold the focal depths of these earthquakes fixed to be unchanged during multi-DD processing by adding in the following constraint equation, in addition to Eq. (1), thusdZi=0(2)where Zi is the focal depth of event i, i = 1, 2, ..., n and n = 8 is the total number of events for which the focal depths have been determined using teleseismic waveform modeling. The joint analysis of local Pg and Sg phases, regional Pn and Sn phases, and teleseismic pP and sP phases with the precise measurements of differential phase arrival times via waveform cross-correlation considerably improves the focal depth determinations.

Seismic tomography

Starting from the relocated earthquakes, we carried out a seismic tomographic inversion to study the three-dimensional Vp and Vs structures of the study region (fig. S5). In addition to the Gorkha earthquake sequence of 533 earthquakes, we also used data from 714 earthquakes that occurred in the study area before the Gorkha earthquake since 1990 (gray circles in Fig. 1). We used regional-scale double-difference tomography (38), which is a generalization of the DD location (35). By taking the curvature of Earth into the Cartesian volume of the grid nodes, we determined the velocity structure for larger scales. The finite-difference method is used for the determination of the travel times and ray paths. These data include 14,893 P-wave and 5831 S-wave arrival times (fig. S1A) recorded at 153 stations (Fig. 1). They form 120,407 P-wave and 22,123 S-wave differential travel times between the events. We used a layered one-dimensional model as an initial model based on velocities from local seismic network (9), regional earthquake relocation and tomography (18, 20), and discontinuities from receiver function analysis (fig. S1B) (3).

Waveform inversions

We constrained absolute focal depths of eight moderate aftershocks (Mw 5.5 to 6.7 in table S1) by modeling the teleseismic waveforms (30°< Δ < 90°) of the direct P and the surface reflection pP and sP waves (35, 36). We used local and regional (0° < Δ < 10°) waveform inversion methods to determine focal mechanisms for 12 smaller earthquakes (Mw < 5.5 in table S2; Fig. 2C) (39). These methods involve matching of complete body waveforms to synthetic waveforms. We assumed that the source can be represented as a point (the centroid) in space and a series of overlapping triangles in time. We calculated synthetic waveforms using layered velocity and density models (fig. S1B). Amplitudes were corrected for geometrical spreading using attenuation t* operator with a value of 1.0 s.


Supplementary material for this article is available at

Fig. S1. Arrival times and the initial P-wave velocity model.

Fig. S2. Bootstrap analysis for evaluating the accuracy of earthquake relocation.

Fig. S3. Comparison between earthquake relocations and the initial locations for the mainshock and 265 selected aftershocks.

Fig. S4. Checkerboard tests of the resulting tomographic images.

Fig. S5. Map view of Vp and Vp/Vs ratio distributions.

Table S1. Locations of 266 earthquakes obtained in this study.

Table S2. Source parameters of 18 earthquakes obtained from earthquake relocation and waveform modeling.

Data file S1. Tomographic solution of Vp and Vp/Vs ratio.

Reference (40)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We used earthquake catalogs from the reviewed International Seismological Centre (ISC) bulletin, China Earthquake Data Center, Nepalese Seismic Network, and waveform data from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center (FDSN code: XQ, 2015-2016) and Data Management Centre of China National Seismic Network at Institute of Geophysics, China Earthquake Administration [SEISDMC; doi:10.11998/SeisDmc/SN (32)]. We are grateful to the teams of China-Nepal border stations, the International Centre for Integrated Mountain Development (ICIMOD), and the NAMASTE project (led by S. Nath Sapkota and J. Nabelek) for conducting the field work. We thank K. Hodges and two anonymous reviewers for their constructive comments, J. Dozal-Cruz and A. Velasco for sharing NAMASTE arrival time picks, and S. Roecker and W. L. Ellsworth for useful discussions. Funding: This research is supported by the grants of Strategic Priority Research Program of Chinese Academy of Sciences (no. XDA20070302), the National Nature Science Foundation of China (nos. 41761144076 and 41490611), the China Scholarship Council, the Cox visiting fellowship of the Stanford University to L.B., the collaborative research program of the Disaster Prevention Research Institute of Kyoto University (no. 29W-03) to L.B. and J.M., and the National Science Foundation (nos. 1545923, 1545933, and 1620646) to M.S.K. and S.L.K., as well as by instrument support from the Kathmandu Center for Research and Education, Chinese Academy of Sciences–Tribhuvan University. Author contributions: L.B. and S.L.K. conceived this study. L.B., S.L.K., and J.M. wrote the manuscript. All authors contributed ideas to the project. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article