Source author record

Mark D. Shattuck

Mark D. Shattuck appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record

Catalog footprint

What is connected

22works
11topics
4close collaborators

Actions

Connect this record

Log in to claim

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

22 published item(s)

preprint2026arXiv

Assessment of scoring functions for computational models of protein-protein interfaces

A goal of computational studies of protein-protein interfaces (PPIs) is to predict the binding site between two monomers that form a heterodimer. The simplest version of this problem is to rigidly re-dock the bound forms of the monomers, which involves generating computational models of the heterodimer and then scoring them to determine the most native-like models. Scoring functions have been assessed previously using rank- and classification-based metrics, however, these methods are sensitive to the number and quality of models in the scoring function training set. We assess the accuracy of seven PPI scoring functions by comparing their scores to a measure of structural similarity to the x-ray crystal structure (i.e. the DockQ score) for a non-redundant set of heterodimers from the Protein Data Bank. For each heterodimer, we generate re-docked models uniformly sampled over DockQ and calculate the Spearman correlation between the PPI scores and DockQ. For some targets, the scores and DockQ are highly correlated; however, for many targets, there are weak correlations. Several physical features can explain the difference between difficult- and easy-to-score targets. For example, strong correlations exist between the score and DockQ for targets with highly intertwined monomers and many interface contacts. We also develop a new score based on only three physical features that matches or exceeds the performance of current PPI scoring functions. These results emphasize that PPI prediction can be improved by focusing on correlations between the PPI score and DockQ and incorporating more discriminating physical features into PPI scoring functions.

preprint2025arXiv

Particle-scale origin of quadrupolar non-affine displacement fields in granular solids

In this work, we identify the local structural defects that control the non-affine displacement fields in jammed disk packings subjected to athermal, quasistatic (AQS) simple shear. While complex non-affine displacement fields typically occur during simple shear, isolated effective quadrupoles are also observed and their probability increases with increasing pressure. We show that the emergence of an isolated effective quadrupole requires the breaking of an interparticle contact that is aligned with low-frequency, spatially extended vibrational modes. Since the Eshelby inhomogeneity problem gives rise to quadrupolar displacement fields in continuum materials, we reformulate and implement Eshelby's equivalent inclusion method (EIM) for jammed disk packings. Using EIM, we show that we can reconstruct the non-affine displacement fields for jammed disk packings in response to applied shear as a sum of discrete Eshelby-like defects that are caused by mismatches in the local stiffnesses of triangles formed from Delaunay triangulation of the disk centers.

preprint2022arXiv

Evolving Programmable Computational Metamaterials

Granular metamaterials are a promising choice for the realization of mechanical computing devices. As preliminary evidence of this, we demonstrate here how to embed Boolean logic gates (AND and XOR) into a granular metamaterial by evolving where particular grains are placed in the material. Our results confirm the existence of gradients of increasing "AND-ness" and "XOR-ness" within the space of possible materials that can be followed by evolutionary search. We measure the computational functionality of a material by probing how it transforms bits encoded as vibrations with zero or non-zero amplitude. We compared the evolution of materials built from mass-contrasting particles and materials built from stiffness-contrasting particles, and found that the latter were more evolvable. We believe this work may pave the way toward evolutionary design of increasingly sophisticated, programmable, and computationally dense metamaterials with certain advantages over more traditional computational substrates.

preprint2022arXiv

The glass-forming ability of binary Lennard-Jones systems

The glass-forming ability (GFA) of alloys, colloidal dispersions, and other particulate materials, as measured by the critical cooling rate $R_c$, can span more than ten orders of magnitude. Even after numerous previous studies, the physical features that control the GFA are still not well understood. For example, it is well-known that mixtures are better glass-formers than monodisperse systems and that particle size and cohesive energy differences among constituents improve the GFA, but it is not currently known how particle size differences couple to cohesive energy differences to determine the GFA. We perform molecular dynamics simulations to determine the GFA of equimolar, binary Lennard-Jones (LJ) mixtures versus the normalized cohesive energy difference $ε_\_$ and mixing energy $\bar ε_{AB}$ between particles A and B. We find several important results. First, the $\log_{10} R_c$ contours in the $\bar ε_{AB}$-$ε_\_$ plane are ellipsoidal for all diameter ratios, and thus $R_c$ is determined by the Mahalanobis distance $d_M$ from a given point in the $\bar ε_{AB}$-$ε_\_$ plane to the center of the ellipsoidal contours. Second, LJ systems for which the larger particles have larger cohesive energy are generally better glass formers than those for which the larger particles have smaller cohesive energy. Third, $d_M(ε_\_,\bar ε_{AB})$ is determined by the relative Voronoi volume difference between particles and local chemical order $S_{AB}$, which gives the average fraction of nearest-neighbor B particles surrounding an A particle and vice-versa. In particular, the shifted Mahalanobis distance $d_M - d^0_M$ versus the shifted chemical order $S_{AB}-S_{AB}^0$ collapses onto a hyperbolic master curve for all diameter ratios.

preprint2020arXiv

Shear response of granular packings compressed above jamming onset

We investigate the mechanical response of jammed packings of repulsive, frictionless spherical particles undergoing isotropic compression. Prior simulations of the soft-particle model, where the repulsive interactions scale as a power-law in the interparticle overlap with exponent $α$, have found that the ensemble-averaged shear modulus $\langle G \rangle$ increases with pressure $P$ as $\sim P^{(α-3/2)/(α-1)}$ at large pressures. However, a deep theoretical understanding of this scaling behavior is lacking. We show that the shear modulus of jammed packings of frictionless, spherical particles has two key contributions: 1) continuous variations as a function of pressure along geometrical families, for which the interparticle contact network does not change, and 2) discontinuous jumps during compression that arise from changes in the contact network. We show that the shear modulus of the first geometrical family for jammed packings can be collapsed onto a master curve: $G^{(1)}/G_0 = (P/P_0)^{(α-2)/(α-1)} - P/P_0$, where $P_0 \sim N^{-2(α-1)}$ is a characteristic pressure that separates the two power-law scaling regions and $G_0 \sim N^{-2(α-3/2)}$. Deviations from this form can occur when there is significant non-affine particle motion near changes in the contact network. We further show that $\langle G (P)\rangle$ is not simply a sum of two power-laws, but $\langle G \rangle \sim (P/P_c)^a$, where $a \approx (α-2)/(α-1)$ in the $P \rightarrow 0$ limit and $\langle G \rangle \sim (P/P_c)^b$, where $b \gtrsim (α-3/2)/(α-1)$ above a characteristic pressure $P_c$. In addition, the magnitudes of both contributions to $\langle G\rangle$ from geometrical families and changes in the contact network remain comparable in the large-system limit for $P >P_c$.

preprint2019arXiv

Intrinsic dissipation mechanisms in metallic glass resonators

Micro- and nano-resonators have important applications including sensing, navigation, and biochemical detection. Their performance is quantified using the quality factor $Q$, which gives the ratio of the energy stored to the energy dissipated per cycle. Metallic glasses are a promising materials class for micro- and nano-scale resonators since they are amorphous and can be fabricated precisely into complex shapes on these lengthscales. To understand the intrinsic dissipation mechanisms that ultimately limit large $Q$-values in metallic glasses, we perform molecular dynamics simulations to model metallic glass resonators subjected to bending vibrations. We calculate the vibrational density of states, redistribution of energy from the fundamental mode of vibration, and $Q$ versus the kinetic energy per atom $K$ of the excitation. In the linear and nonlinear response regimes where there are no atomic rearrangements, we find that $Q \rightarrow \infty$ (since we do not consider coupling to the environment). We identify a characteristic $K_r$ above which atomic rearrangements occur, and there is significant energy leakage from the fundamental mode to higher frequencies, causing finite $Q$. Thus, $K_r$ is a critical parameter determining resonator performance. We show that $K_r$ decreases as a power-law, $K_r\sim N^{-k},$ with increasing system size $N$, where $k \approx 1.3$. We estimate the critical strain $\langle γ_r \rangle \sim 10^{-8}$ for micron-sized resonators below which atomic rearrangements do not occur, and thus large $Q$-values can be obtained when they are operated below $γ_r$. We find that $K_r$ for amorphous resonators is comparable to that for resonators with crystalline order.

preprint2019arXiv

Pressure-dependent shear response of jammed packings of spherical particles

The mechanical response of packings of purely repulsive, spherical particles to athermal, quasistatic simple shear near jamming onset is highly nonlinear. Previous studies have shown that, at small pressure $p$, the ensemble-averaged static shear modulus $\langle G-G_0 \rangle$ scales with $p^α$, where $α\approx 1$, but above a characteristic pressure $p^{**}$, $\langle G-G_0 \rangle \sim p^β$, where $β\approx 0.5$. However, we find that the shear modulus $G^i$ for an individual packing typically decreases linearly with $p$ along a geometrical family where the contact network does not change. We resolve this discrepancy by showing that, while the shear modulus does decrease linearly within geometrical families, $\langle G \rangle$ also depends on a contribution from discontinuous jumps in $\langle G \rangle$ that occur at the transitions between geometrical families. For $p > p^{**}$, geometrical-family and rearrangement contributions to $\langle G \rangle$ are of opposite signs and remain comparable for all system sizes. $\langle G \rangle$ can be described by a scaling function that smoothly transitions between the two power-law exponents $α$ and $β$. We also demonstrate the phenomenon of {\it compression unjamming}, where a jammed packing can unjam via isotropic compression.

preprint2016arXiv

Reconstruction of Ordinary Differential Equations From Time Series Data

We develop a numerical method to reconstruct systems of ordinary differential equations (ODEs) from time series data without {\it a priori} knowledge of the underlying ODEs using sparse basis learning and sparse function reconstruction. We show that employing sparse representations provides more accurate ODE reconstruction compared to least-squares reconstruction techniques for a given amount of time series data. We test and validate the ODE reconstruction method on known 1D, 2D, and 3D systems of ODEs. The 1D system possesses two stable fixed points; the 2D system possesses an oscillatory fixed point with closed orbits; and the 3D system displays chaotic dynamics on a strange attractor. We determine the amount of data required to achieve an error in the reconstructed functions to less than $0.1\%$. For the reconstructed 1D and 2D systems, we are able to match the trajectories from the original ODEs even at long times. For the 3D system with chaotic dynamics, as expected, the trajectories from the original and reconstructed systems do not match at long times, but the reconstructed and original models possess similar Lyapunov exponents. Now that we have validated this ODE reconstruction method on known models, it can be employed in future studies to identify new systems of ODEs using time series data from deterministic systems for which there is no currently known ODE model.

preprint2015arXiv

Asymmetric crystallization during cooling and heating in model glass-forming systems

We perform molecular dynamics (MD) simulations of the crystallization process in binary Lennard-Jones systems during heating and cooling to investigate atomic-scale crystallization kinetics in glass-forming materials. For the cooling protocol, we prepared equilibrated liquids above the liquidus temperature $T_l$ and cooled each sample to zero temperature at rate $R_c$. For the heating protocol, we first cooled equilibrated liquids to zero temperature at rate $R_p$ and then heated the samples to temperature $T > T_l$ at rate $R_h$. We measured the critical heating and cooling rates $R_h^*$ and $R_c^*$, below which the systems begin to form a substantial fraction of crystalline clusters during the heating and cooling protocols. We show that $R_h^* > R_c^*$, and that the asymmetry ratio $R_h^*/R_c^*$ includes an intrinsic contribution that increases with the glass-forming ability (GFA) of the system and a preparation-rate dependent contribution that increases strongly as $R_p \rightarrow R_c^*$ from above. We also show that the predictions from classical nucleation theory (CNT) can qualitatively describe the dependence of the asymmetry ratio on the GFA and preparation rate $R_p$ from the MD simulations and results for the asymmetry ratio measured in Zr- and Au-based bulk metallic glasses (BMG). This work emphasizes the need for and benefits of an improved understanding of crystallization processes in BMGs and other glass-forming systems.

preprint2015arXiv

Beyond packing of hard spheres: The effects of core softness, non-additivity, intermediate-range repulsion, and many-body interactions on the glass-forming ability of bulk metallic glasses

When a liquid is cooled well below its melting temperature at a rate that exceeds the critical cooling rate $R_c$, the crystalline state is bypassed and an amorphous glassy state forms instead. $R_c$ (or the corresponding critical casting thickness $d_c$) characterizes the glass-forming ability (GFA) of each material. While silica is an excellent glass-former with small $R_c<10^{-2}$ K/s, pure metals and most alloys are poor glass-formers with large $R_c>10^{10}$ K/s. Only in the past thirty years have bulk metallic glasses (BMGs) been identified with $R_c$ approaching that for silica. Recent simulations have shown that hard-sphere models are able to identify the atomic size ratio and number fraction regime where BMGs exist with critical cooling rates more than 13 orders of magnitude smaller than those for pure metals. However, there are many other features of interatomic potentials beyond hard-core interactions. How do these other features affect the glass-forming ability of BMGs? We perform molecular dynamics simulations to determine how variations in the softness and non-additivity of the repulsive core and form of the interatomic pair potential at intermediate distances affect the GFA of binary alloys. These variations in the interatomic pair potential allow us to introduce geometric frustration and change the crystal phases that compete with glass formation. We also investigate the effect of tuning the strength of the many-body interactions from zero to the full embedded atom model on the GFA for pure metals. We then employ the full embedded atom model for binary BMGs and show that hard-core interactions play the dominant role in setting the GFA of alloys, while other features of the interatomic potential only change the GFA by one to two orders of magnitude.

preprint2015arXiv

Onset and cessation of motion in hydrodynamically sheared granular beds

We performed molecular dynamics simulations of granular beds driven by a model hydrodynamic shear flow to elucidate general grain-scale mechanisms that determine the onset and cessation of sediment transport. By varying the Shields number (the nondimensional shear stress at the top of the bed) and particle Reynolds number (the ratio of particle inertia to viscous damping), we explore how variations of the fluid flow rate, particle inertia, and fluid viscosity affect the onset and cessation of bed motion. For low to moderate particle Reynolds numbers, a critical boundary separates mobile and static states. Transition times between these states diverge as this boundary is approached both from above and below. At high particle Reynolds number, inertial effects become dominant, and particle motion can be sustained well below flow rates at which mobilization of a static bed occurs. We also find that the onset of bed motion (for both low and high particle Reynolds numbers) is described by Weibullian weakest-link statistics, and thus is crucially dependent on the packing structure of the granular bed, even deep beneath the surface.

preprint2015arXiv

Outcome prediction in mathematical models of immune response to infection

Clinicians need to predict patient outcomes with high accuracy as early as possible after disease inception. In this manuscript, we show that patient-to-patient variability sets a fundamental limit on outcome prediction accuracy for a general class of mathematical models for the immune response to infection. However, accuracy can be increased at the expense of delayed prognosis. We investigate several systems of ordinary differential equations (ODEs) that model the host immune response to a pathogen load. Advantages of systems of ODEs for investigating the immune response to infection include the ability to collect data on large numbers of `virtual patients', each with a given set of model parameters, and obtain many time points during the course of the infection. We implement patient-to-patient variability $v$ in the ODE models by randomly selecting the model parameters from Gaussian distributions with variance $v$ that are centered on physiological values. We use logistic regression with one-versus-all classification to predict the discrete steady-state outcomes of the system. We find that the prediction algorithm achieves near $100\%$ accuracy for $v=0$, and the accuracy decreases with increasing $v$ for all ODE models studied. The fact that multiple steady-state outcomes can be obtained for a given initial condition, i.e. the basins of attraction overlap in the space of initial conditions, limits the prediction accuracy for $v>0$. Increasing the elapsed time of the variables used to train and test the classifier, increases the prediction accuracy, while adding explicit external noise to the ODE models decreases the prediction accuracy. Our results quantify the competition between early prognosis and high prediction accuracy that is frequently encountered by clinicians.

preprint2014arXiv

Can the packing efficiency of binary hard spheres explain the glass-forming ability of bulk metallic glasses?

We perform molecular dynamics simulations to compress binary hard spheres into jammed packings as a function of the compression rate $R$, size ratio $α$, and number fraction $x_S$ of small particles to determine the connection between the glass-forming ability (GFA) and packing efficiency in bulk metallic glasses (BMGs). We define the GFA by measuring the critical compression rate $R_c$, below which jammed hard-sphere packings begin to form "random crystal" structures with defects. We find that for systems with $α\gtrsim 0.8$ that do not de-mix, $R_c$ decreases strongly with $Δϕ_J$, as $R_c \sim \exp(-1/Δϕ_J^2)$, where $Δϕ_J$ is the difference between the average packing fraction of the amorphous packings and random crystal structures at $R_c$. Systems with $α\lesssim 0.8$ partially de-mix, which promotes crystallization, but we still find a strong correlation between $R_c$ and $Δϕ_J$. We show that known metal-metal BMGs occur in the regions of the $α$ and $x_S$ parameter space with the lowest values of $R_c$ for binary hard spheres. Our results emphasize that maximizing GFA in binary systems involves two competing effects: minimizing $α$ to increase packing efficiency, while maximizing $α$ to prevent de-mixing.

preprint2014arXiv

The glass-forming ability of model metal-metalloid alloys

Bulk metallic glasses (BMGs) are amorphous alloys with desirable mechanical properties and processing capabilities. To date, the design of new BMGs has largely employed empirical rules and trial-and-error experimental approaches. Ab initio computational methods are currently prohibitively slow to be practically used in searching the vast space of possible atomic combinations for bulk glass formers. Here, we perform molecular dynamics simulations of a coarse-grained, anisotropic potential, which mimics interatomic covalent bonding, to measure the critical cooling rates for metal-metalloid alloys as a function of the atomic size ratio $σ_S/σ_L$ and number fraction $x_S$ of the metalloid species. We show that the regime in the space of $σ_S/σ_L$ and $x_S$ where well-mixed, optimal glass formers occur for patchy and LJ particle mixtures coincides with that for experimentally observed metal-metalloid glass formers. Our simple computational model provides the capability to perform combinatorial searches to identify novel glass-forming alloys.

preprint2014arXiv

The statistics of frictional families

We develop a theoretical description for mechanically stable frictional packings in terms of the difference between the total number of contacts required for isostatic packings of frictionless disks and the number of contacts in frictional packings, $m=N_c^0-N_c$. The saddle order $m$ represents the number of unconstrained degrees of freedom that a static packing would possess if friction were removed. Using a novel numerical method that allows us to enumerate disk packings for each $m$, we show that the probability to obtain a packing with saddle order $m$ at a given static friction coefficient $μ$, $P_m(μ)$, can be expressed as a power-series in $μ$. Using this form for $P_m(μ)$, we quantitatively describe the dependence of the average contact number on friction coefficient for static disk packings obtained from direct simulations of the Cundall-Strack model for all $μ$ and $N$.

preprint2013arXiv

Computational studies of the glass-forming ability of model bulk metallic glasses

Bulk metallic glasses (BMGs) are produced by rapidly thermally quenching supercooled liquid metal alloys below the glass transition temperature at rates much faster than the critical cooling rate R_c below which crystallization occurs. The glass-forming ability of BMGs increases with decreasing R_c, and thus good glass-formers possess small values of R_c. We perform molecular dynamics simulations of binary Lennard-Jones (LJ) mixtures to quantify how key parameters, such as the stoichiometry, particle size difference, attraction strength, and heat of mixing, influence the glass-formability of model BMGs. For binary LJ mixtures, we find that the best glass-forming mixtures possess atomic size ratios (small to large) less than 0.92 and stoichiometries near 50:50 by number. In addition, weaker attractive interactions between the smaller atoms facilitate glass formation, whereas negative heats of mixing (in the experimentally relevant regime) do not change R_c significantly. These studies represent a first step in the development of computational methods for quantitatively predicting glass-formability.

preprint2013arXiv

Isostaticity at Frictional Jamming

Amorphous packings of frictionless, spherical particles are isostatic at jamming onset, with the number of constraints (contacts) equal to the number of degrees of freedom. Their structural and mechanical properties are controlled by the interparticle contact network. In contrast, amorphous packings of frictional particles are typically hyperstatic at jamming onset. We perform extensive numerical simulations in two dimensions of the geometrical asperity (GA) model for static friction, to further investigate the role of isostaticity. In the GA model, interparticle forces are obtained by summing up purely repulsive central forces between periodically spaced circular asperities on contacting grains. We compare the packing fraction, contact number, mobilization distribution, and vibrational density of states using the GA model to those generated using the Cundall-Strack (CS) approach. We find that static packings of frictional disks obtained from the GA model are mechanically stable and isostatic when we consider interactions between asperities on contacting particles. The crossover in the structural and mechanical properties of static packings from frictionless to frictional behavior as a function of the static friction coefficient coincides with a change in the type of interparticle contacts and the disappearance of a peak in the density of vibrational modes for the GA model. These results emphasize that mesoscale features of the model for static friction play an important role in determining the properties of granular packings.

preprint2013arXiv

Particle-scale reversibility in athermal particulate media below jamming

We perform numerical simulations of athermal repulsive frictionless disks and spheres in two and three spatial dimensions undergoing cyclic quasi-static simple shear to investigate particle-scale reversible motion. We identify three classes of steady-state dynamics as a function of packing fraction ϕand maximum strain amplitude per cycle γ_{\rm max}. Point-reversible states, where particles do not collide and exactly retrace their intra-cycle trajectories, occur at low ϕand γ_{\rm max}. Particles in loop-reversible states undergo numerous collisions and execute complex trajectories, but return to their initial positions at the end of each cycle. Loop-reversible dynamics represents a novel form of self-organization that enables reliable preparation of configurations with specified structural and mechanical properties over a broad range of ϕfrom contact percolation to jamming onset at ϕ_J. For sufficiently large ϕand γ_{\rm max}, systems display irreversible dynamics with nonzero self-diffusion.

preprint2013arXiv

Vibrations in jammed solids: Beyond linear response

We propose a `phase diagram' for particulate systems that interact via purely repulsive contact forces, such as granular media and colloidal suspensions. We identify and characterize two distinct classes of behavior as a function of the input kinetic energy per degree of freedom $T_0$ and packing fraction deviation above and below jamming onset $Δϕ=ϕ- ϕ_J$ using numerical simulations of purely repulsive frictionless disks. Iso-coordinated solids (ICS) only occur above jamming for $Δϕ> Δϕ_c(T_0)$; they possess average coordination number equal to the isostatic value ($< z> = z_{\rm iso}$) required for mechanically stable packings. ICS display harmonic vibrational response, where the density of vibrational modes from the Fourier transform of the velocity autocorrelation function is a set of sharp peaks at eigenfrequencies $ω_k^d$ of the dynamical matrix evaluated at $T_0=0$. Hypo-coordinated solids (HCS) occur both above and below jamming onset within the region defined by $Δϕ> Δϕ^*_-(T_0)$, $Δϕ< Δϕ^*_+(T_0)$, and $Δϕ> Δϕ_{cb}(T_0)$. In this region, the network of interparticle contacts fluctuates with $< z> \approx z_{\rm iso}/2$, but cage-breaking particle rearrangements do not occur. The HCS vibrational response is nonharmonic, {\it i.e} the density of vibrational modes $D(ω)$ is not a collection of sharp peaks at $ω_k^d$, and its precise form depends on the measurement method. For $Δϕ> Δϕ_{cb}(T_0)$ and $Δϕ< Δϕ^*_{-}(T_0)$, the system behaves as a hard-particle liquid.

preprint2013arXiv

Vibrations of Jammed Disk Packings with Hertzian Interactions

Contact breaking and Hertzian interactions between grains can both give rise to nonlinear vibrational response of static granular packings. We perform molecular dynamics simulations at constant energy in 2D of frictionless bidisperse disks that interact via Hertzian spring potentials as a function of energy and measure directly the vibrational response from the Fourier transform of the velocity autocorrelation function. We compare the measured vibrational response of static packings near jamming onset to that obtained from the eigenvalues of the dynamical matrix to determine the temperature above which the linear response breaks down. We compare packings that interact via single-sided (purely repulsive) and double-sided Hertzian spring interactions to disentangle the effects of the shape of the potential from contact breaking. Our studies show that while Hertzian interactions lead to weak nonlinearities in the vibrational behavior (e.g. the generation of harmonics of the eigenfrequencies of the dynamical matrix), the vibrational response of static packings with Hertzian contact interactions is dominated by contact breaking as found for systems with repulsive linear spring interactions.

preprint2011arXiv

The contact percolation transition

Typical quasistatic compression algorithms for generating jammed packings of athermal, purely repulsive particles begin with dilute configurations and then apply successive compressions with relaxation of the elastic energy allowed between each compression step. It is well-known that during isotropic compression athermal systems with purely repulsive interactions undergo a jamming transition at packing fraction $ϕ_J$ from an unjammed state with zero pressure to a jammed, rigid state with nonzero pressure. Using extensive computer simulations, we show that a novel second-order-like transition, the contact percolation transition, which signals the formation of a system-spanning cluster of mutually contacting particles, occurs at $ϕ_P < ϕ_J$, preceding the jamming transition. By measuring the number of non-floppy modes of the dynamical matrix, and the displacement field and time-dependent pressure following compression, we find that the contact percolation transition also heralds the onset of complex spatiotemporal response to applied stress. Thus, highly heterogeneous, cooperative, and non-affine particle motion occurs in unjammed systems significantly below the jamming transition for $ϕ_P < ϕ< ϕ_J$, not only for jammed systems with $ϕ> ϕ_J$.