Deciphering and steering population-level response under spatial drug heterogeneity on microhabitat structures

Deciphering and steering population-level response under spatial drug heterogeneity on microhabitat structures | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Deciphering and steering population-level response under spatial drug heterogeneity on microhabitat structures View ORCID Profile Zhijian Hu , View ORCID Profile Kevin Wood doi: https://doi.org/10.1101/2025.02.13.638200 Zhijian Hu 1 Department of Biophysics, University of Michigan , Ann Arbor, USA 2 Department of Mathematics, University of Michigan , Ann Arbor, USA 3 Center for the Study of Complex Systems, University of Michigan , Ann Arbor, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Zhijian Hu For correspondence: zhijianh{at}umich.edu Kevin Wood 1 Department of Biophysics, University of Michigan , Ann Arbor, USA 3 Center for the Study of Complex Systems, University of Michigan , Ann Arbor, USA 4 Department of Physics, University of Michigan , Ann Arbor, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Kevin Wood Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Bacteria and cancer cells inhabit spatially heterogeneous environments, where migration shapes microhabitat structures critical for colonization and metastasis. The interplay between growth, migration, and spatial structure complicates the prediction of population responses to drug treatment—such as clearance or persistence—even under the same spatially averaged growth rate. Accurately predicting these responses is essential for designing effective treatment strategies. Here, we propose a minimal growth-migration model to study population dynamics on discrete microhabitat structures under spatial drug heterogeneity. By applying a kernel transformation, we map the original structure to an effective fully connected graph and derive a new exact criterion for population response based on a regularized Laplacian kernel reweighted by local growth rates. This criterion connects to forest closeness centrality and yields analytical bounds and sufficient conditions for population growth or decline. We find that higher structural connectivity—like increased migration—generally promotes decline. Our framework also informs optimal spatial drug assignments, which reduce to selecting interconnected subcores in the effective complete graph. For partially controllable microhabitats or unknown drug distributions, we identify strategies that ensure population decline. Overall, our results offer a new theoretical perspective on drug response in spatially structured populations and provide practical guidance for optimizing spatially explicit dosing strategies in heterogeneous environments. Introduction Healthcare-related cell communities, such as bacterial infections and cancer, have increasingly been recognized as complex ecosystems over the past decades [ 1 – 3 ]. These systems exhibit emergent collective behaviors, including antibiotic resistance and cancer invasion, which arise from the interplay of ecological and evolutionary processes. The complex interactions between different phenotypes and genotypes have inspired novel therapeutic approaches, such as containment strategies and adaptive therapies [ 4 – 8 ], which leverage competition between cell types. From a dynamical systems perspective, drugs constitute an integral part of the environment, and their interactions with cells—particularly under spatially heterogeneous drug distributions—have gained attention in recent years [ 5 , 9 – 39 ]. At the evolutionary timescale, studies on cancer have demonstrated that spatial drug heterogeneity can accelerate the evolution of drug resistance [ 10 ], shortening the time between initial treatment and the emergence of resistant clones at detectable levels. This phenomenon has been captured in minimal two-well models and extended ladder models, where differences in drug concentration modulate growth rates [ 16 , 17 ]. Similar acceleration effects have been observed in bacterial populations through laboratory experiments [ 11 , 12 ], though recent theoretical work using mean first passage time models suggests that spatial heterogeneity can either accelerate or decelerate resistance evolution depending on mutation and migration rates [ 10 , 36 , 39 , 40 ]. These contrasting outcomes underscore the complexity of resistance evolution and suggest that optimal spatial drug heterogeneity design could help slow resistance development. At the ecological timescale, similar variability is observed in patient outcomes. Even under identical drug dosages, patients exhibit different levels of pathogen or cancer cell clearance between the first treatment and the final day of therapy. For instance, despite a standard 14-day antibiotic treatment regimen, approximately 20% of patients fail to achieve sufficient clearance of H. pylori [ 41 ]. Similarly, cancer patients face recurrence risks due to insufficient clearance [ 5 , 42 , 43 ], which may be influenced by spatial heterogeneities in drug absorption and pharmacokinetics/pharmacodynamics (PK/PD) [ 25 ]. Unlike well-mixed laboratory conditions, the in vivo tumor and microbial environments exhibit spatial complexity, making traditional well-mixed drug-dose response measurements insufficient for predicting clinical outcomes. Recent work [ 44 ] has modeled population-level responses as an eigenvalue problem, where cell survival or extinction is governed by the interplay of spatial growth, migration, and microhabitat connectivity in a one-dimensional setting. This framework has been extended to model the emergence of multidrug resistance under natural selection [ 45 ]. However, the entanglement of structural, migratory, and growth-related factors in these eigenvalue formulations limits their interpretability. Additionally, clinical applications require extending these results beyond one-dimensional models to account for multiple discrete microhabitat scenarios such as cancer metastasis [ 46 – 51 ], lymph node networks [ 52 – 54 ], and bacterial translocation across body compartments [ 15 , 55 , 56 ]. To better understand how spatial drug heterogeneity shapes population-level responses—specifically, whether cell populations grow or decline—we develop a minimal growth-migration model defined on discrete microhabitat networks. In this model, local growth rates are constrained by drug-induced limits, ranging from a maximum drug-free growth rate to a minimum death rate. Migration allows for population movement across interconnected microhabitats, embedding spatial structure into the dynamics. We theoretically analyze this model by transforming the original network into an effective fully connected graph using a regularized Laplacian kernel. This transformation enables the derivation of an exact criterion for determining population response, reweighted by local growth rates and interpreted through a centrality measure known as forest closeness centrality. This centrality interpretation provides an efficient and analytically tractable condition for predicting whether a population will persist or be cleared under a given spatial drug distribution. Interestingly, we find that increasing global structural diffusion, such as through enhanced migration or more inter-microhabitat connectivities, leads to a smooth transition from population growth to decline when the spatial drug heterogeneity is fixed. For a given microhabitat structure, population dynamics are not only governed by the spatially averaged growth rate and the migration rate, but the spatial arrangement of growth rates still plays a critical role. Our framework informs optimal spatial drug assignments, which reduce to selecting interconnected subcores in the effective complete graph. For partially controllable microhabitats, we can design optimized strategies for population clearance; for unknown spatial drug heterogeneities, under a centrality-based heuristic, we identify parameters that ensure robust population decline, independent of spatial drug arrangements. These findings have important implications for understanding and manipulating drug responses in complex biological systems. Health-related applications—such as bacterial infections and metastatic cancer—often involve heterogeneous spatial environments where growth, migration, and structure interact in nontrivial ways. Our framework offers a generalizable and analytically grounded approach to assess when populations will decline in such settings, providing a new lens through which to interpret treatment outcomes. Moreover, the model serves as a foundation for future extensions, including those incorporating temporal fluctuations, environmental feedback, or interspecies interactions. Clinically, this work highlights the potential of spatially explicit dosing strategies, tailored not only to drug strength but also to the structure of the tissue or infection site, to improve treatment outcomes for both infections and solid tumors. Basic set-up of model system Cell populations like bacteria, or cancer cells, can actively or passively move between different sites by migration or metastasis. When accounting for fluctuations in population growth and the stochastic nature of drug-induced clearance [ 57 ], the system can be described using stochastic microscopic dynamics governed by a master equation (see SI). However, in human bodies, where bacterial and cancer cell populations typically exceed thousands, we can apply the Van Kampen approximation [ 58 ]. This allows us to ignore microscopic dynamic details and derive a simplified macroscopic growth-migration equation This is equivalent to a linearized exponential growth from logistic growth, by ignoring the competition or density-dependent effect, where is a diagonal matrix describing growth rates induced by drug at different microhabitats. And L is the Laplacian matrix describing the multi-microhabitat structure. β describes the migration rate and u = ( u 1 , u 2 , …, u N ) T describes the population densities on different microhabitats. N is the total number of microhabitats. Since drug-induced growth rates { g i } cannot be infinite, we require g i ∈ [ d 0 , g 0 ], where d 0 0 is the maximum growth rate or drug-free growth rate, depending on drug type and cell type, or other environmental factors like nutrition [ 59 – 61 ]. For example, single bactericidal drug like Ampicillin can be fitted into a hill fuction to get maximum growth rate at a drug-free environment, or the maximum death rate when the AMP concentration is extremely high [ 62 ]. Equation (1) is the minimal growth-migration model required to understand short-term population response with different spatial drug heterogeneities and different multi-microhabitat strcutures. Since the drug concentration must be sufficiently high to eliminate cancer or bacterial cells, we require that the spatially averaged growth rate satisfies . For population decline or clearance, the condition Ω = G − βL ≺ 0 must hold, or equivalently, the largest eigenvalue λ 0 = ∥Ω∥ = ∥ G − βL ∥ < 0. However, solving this exactly is typically intractable. Applying first-order perturbation theory, as in recent studies [ 44 , 45 ], we approximate the largest eigenvalue as Since , we obtain ⟨ u 0 | δ g | u 0 ⟩ = 0, indicating that this approximation only captures well-mixed effects while neglecting spatial drug distribution and multi-microhabitat structure. Higher-order approximations may be necessary but are challenging to compute. Results Kernel transformation and new criterion for population response Instead of relying on perturbation theory, we derive an exact solution for the clearance or decline criterion based on a “kernel” transformation. For the kernel transformation, we first decompose G as , where I is the identity matrix and is an orthogonal diagonal matrix with the ith diagonal entry as , representing the relative growth information over microhabitats compared to the maximum death rate | d 0 |. Based on Schur complements (see SI for details), we show that the condition Ω = G − βL ≺ 0 is equivalent to where K = ( I + γL ) −1 is the regularized Laplacian kernel matrix naturally emerging from our minimal growth-migration dynamics, and is the ratio of migration rate to maximum death rate. Traditionally, the parameter γ characterizes the extent of diffusion or spread over the graph structure. Here, it effectively acts as a rescaled migration intensity, modulating how spatial structure shapes population outcomes. The matrix K , also known in some contexts as the forest matrix or the parametrized matrix forest index , has appeared in a wide range of network applications. It is at least row-stochastic (and becomes doubly stochastic for undirected graphs), mapping the original sparse structure into an effective, fully connected graph. Off-diagonal entries K ij have been used to quantify similarity between spatial locations in disease-spread models, such as COVID-19 transmission networks [ 63 ]. More broadly, K has been applied for link prediction in diverse domains, including social and biological networks [ 64 ], protein–protein interaction (PPI) networks [ 65 , 66 ], and virus–host interaction inference [ 67 ]. A particularly relevant measure derived from K is the forest closeness centrality [ 68 , 69 ], defined as which approximately corresponds to the inverse of the diagonal element K ii . For simplicity, we will refer to as the centrality measure in this paper. While this matrix has traditionally been engineered for tasks such as node ranking, link prediction, or similarity assessment, our framework reveals its natural appearance as a kernel governing the interplay between growth and migration. This connection provides a simple but novel interpretation of centrality as a predictor of population response and optimized clearance strategies, which we further illuminate in the subsequent sections. Rewrite where ∥ · ∥ is chosen as the infinite norm | · | ∞ , the maximum row sum of absolute value of matrix A . here can be interpreted as the length-k walk from microhabitat i to microhabitat j . So now λ max ( A ) becomes a special norm measuring the maximum total walk length from a microhabitat i to every microhabitat including i itself (see SI for details). Again, this tells us the new largest eigenvalue is capturing global information of the original multi-microhabitat structure. It’s hard to find out the analytical solution of λ max ( A ) due to the inhomgeneous nature of growth-weight distance between microhabitats. However, based on this walk-length interpretation, we find out the lower and upper bound for λ max ( A ), where K ii 1 /C i in the lower bound represents the inverse of forest closeness centrality, and the upper bound is the maximum row sum of length-1 walk (see SI for details). We can also treat each as the edge value of the efficient fully connected graph.Then these 2 bounds inform that, finding λ max ( A ) here is approximately equivalent to identifying key microhabitats with maximum growth-weighted centrality, or with maximum weighted connectivities. λ max ( A ) = max i K ii δ i holds when microhabitat i becomes effectively isolated; holds when the structure has local symmetry and A becomes a circulant matrix(see SI for details). Either K ii or K ij contains global information of the original structure due to the kernel transformation. If we do the same calculation for Ω = G − βL instead of A , we can only find out g i − βk ( i ) ≤ λ max (Ω) ≤ g i , where k ( i ) is the degree of ith mcriohabitat. These are looser bounds only containing localized information like degree. So by doing the kernel formation, we can find out a tighter bound of λ max with meaningful interpretations. If or max i K ii δ i > 1, we thus find out sufficient conditions for population decline or population growth (or , max i K ii δ i < 1 are necessary conditions for population growth, population decline). These 2 criteria can be applied to quickly determine population response without directly calculating λ max ( A ) itself. While is in a first-order form considering interaction effect, max i K ii δ i has zero-order form considering only centrality effect, decoupling growth information and structure information 1 . So in practice, we can always perform a quick diagnostic check of the population response using the quantity max i K ii δ i . If max i K ii δ i > 1, we can immediately predict population growth. Conversely, if max i K ii δ i < 1, this provides a necessary condition for population decline. As a second step, we evaluate the upper bound ; if this quantity is also smaller than 1, we can conclude that population decline will occur. Since the diagonal element K ii , also known as the forest closeness centrality C i , serves as a useful and interpretable indicator under growth-migration dynamics, we propose to refer to K ii as the inverse dynamic-related centrality , and denote its reciprocal C i as the dynamic-related centrality . In this study, we primarily focus on how inverse centrality values K ii can be leveraged to predict population responses. Because these responses are influenced by two spatial effects — spatial drug heterogeneity (i.e., the spatial growth distribution { δ i }) and the spatial structure of the system — we organize our analysis as follows: Given a fixed spatial drug heterogeneity { δ i }, we analyze how different microhabitat structures, as captured by the matrix K , influence population-level responses, and examine the role that K ii plays in shaping these dynamics. Given a known multi-microhabitat structure (i.e., a fixed matrix K ), we explore how different growth distributions { δ i } affect the outcome, and how this insight can be used to optimize treatment strategies by exploiting structural information — particularly the inverse centralities K ii . Microhabitat structure effects under given spatial drug heterogeneity We aim to investigate how microhabitat structures influence population responses under a fixed spatial drug heterogeneity. Specifically, we address the following questions: 1. How do different microhabitat structures with the same connectivity affect population responses? 2. How does varying the level of connectivity alter spatial drug responses across different structural configurations? Overall, our goal is to understand how structural effects on population dynamics can be predicted by the dynamic-related centrality C i , or equivalently, by its inverse K ii . Distinct population responses induced by different multi-microhabitat structures predicted by kernel criteria Starting from a fixed spatial drug heterogeneity (see Figure 2A ), we demonstrate how different multi-microhabitat structures lead to distinct population responses, which can be precisely captured by checking max i K ii δ i and . The inverse value of dynamic-related centralities C are visualized for two different multi-microhabitat structures, with color gradients indicating different values (see Figure 1B ). All values of inversed C (or K ii ) are upper bounded by 1 (see SI). Different transparencies represent varying growth rate values. Download figure Open in new tab Figure 1. Illustration of Population Response Under Drug Exposure. Example of multiple bacterial microhabitats and their diverging population responses under high bactericidal drug doses. Despite receiving the same drug dose, bacterial populations may exhibit microhabitat-dependent responses, even when the population would decline in a well-mixed environment. Download figure Open in new tab Figure 2. Different microhabitat structures with the same spatial growth distributions induce different population responses, explained by our criterion. A . Example distribution of spatial growth rates across 10 connected microhabitats. Growth rates are bounded by the drug-related maximum growth rate and death rate, represented by grey dashed lines. B . Two different microhabitat structures with the same number of microhabitats and migration connections, along with their corresponding K ii δ i values shown as bar plots. In each structure, colors indicate K ii , the inverse of dynamic-related centrality values C i , while different transparencies represent relative growth δ i . The microhabitat with the highest K ii δ i in each structure is highlighted with a triangle, colored light yellow for population growth and light blue for population decline, as confirmed by the adjacent bar plots. Notably, the microhabitat with the highest K ii does not necessarily correspond to the highest K ii δ i . The yellow check mark near the top bar plot is showing that the population growth can be predicted by checking max i K ii δ i > 1. The light blue cross mark near the bottom bar plot indicates that max i K ii δ i < 1 serves as a necessary condition for population declinew, and further check of is required. C . Two spatial dynamic configurations demonstrating different population responses. The grey dashed line represents the initial population density. As time increases, the density curves change from grey to dark blue. The dynamics induced by structure 1 exhibit a growth trend, consistent with our theoretical prediction, while those induced by structure 2 show a decline trend predicted by . For both 2 structures, β = 1, g 0 = 2, d 0 = −5, ⟨ g ⟩ = −1.5, and the edge number is fixed at | E | = 10. Notably, the locations of g max (or δ max ), max K ii , and max i K ii δ i do not necessarily coincide. We can approximate K ii δ i directly from the visualized networks by observing the color and transparency of microhabitats(nodes), which are highlighted in different colors. The adjacent bar plots reveal that in one structure, at least one K ii δ i value exceeds 1, whereas in the other, all values remain below 1. Consequently, the first multi-microhabitat structure leads to population growth, while the second indicatively results in decline. It’s proved by checking . The population decline is as illustrated by the spatial-temporal dynamics in Figure 2C . Our simple criterion thus serves as a powerful predictive tool for population responses without requiring detailed dynamic simulations. Interestingly, dynamic-related centrality C itself, is proportional to node degree, correlating with many conventional centrality measures (see SI), although the importance rankings are not exactly the same. It bears strong similarity to closeness centrality. In the population responses of 2 example structures we haven shown, structure 2 with more even connectivities for different microhabitats, tend to decline compared to structure 1. This can be explained from a centrality perspective: high C usually represents good connectivities with other microhabitats, if under a relatively high-drug dose regime as shown in Figure 2A , the high C tends to connect with more “sink” microhabitats with negative growth rates, increasing the probability for population decline. In structure 1, the microhabitat 4, with positive growth rate high inversed centrality value (low centrality C value), is pretty isolated with only one connection to other microhabitats, thus increasing the chance of survival. Increased connectivity accelerates centrality-predicted population decline Since dynamic-related centrality reflects the importance of edge number, and the second structure with a more even connectivity in Figure 2 leads to population decline, it suggests that increasing global diffusion over different microhabitats can effectively decrease K ii , thus mitigating population growth by max K ii δ i < 1. Increasing migration rate can be one way to enhance diffusion(we can also see that from the mathmatical form of , whee K ii ∝ β −1 ). Increasing the average number of edges or connection density may be another way to enhance global diffusion thu increasing the probability of population decline. To investigate this, we tune the edge density , defined as the ratio of the total number of edges in a given graph to the number of edges in a complete graph, and examine its effect on max K ii δ i . Using 1000 stochastic samples from random graph for each edge density, we observe that as edge density increases, sample-averaged max K ii δ i decreases rapidly and eventually crosses the critical threshold max K ii δ i = 1 (see Figure 3 ). Figure 3A illustrates four different multi-microhabitat structures with increasing edge densities, where the colors gradually fade to white, indicating decreasing max K ii δ i values. This finding suggests that densely connected microhabitats may require lower drug doses to induce population decline, whereas sparsely connected structures may necessitate higher drug doses to achieve the same effect. Download figure Open in new tab Figure 3. Increasing edge/connection density enhances the tendency for population decline. A. Four example microhabitat structures(microhabitat number N = 40) with varying connection densities but the same spatial growth rate distribution. Color gradients represent dynamic-related centrality values C , while different transparencies indicate relative growth values δ . As connectivity increases, the colors gradually fade to white, indicating lower centrality values. B . The maximum Cδ value decreases as connection density ρ increases. As the curve crosses the critical threshold max Cδ = 1, population response shifts from growth to decline. Connectivity has a similar effect to migration rate (see SI). The curve represents an average over 1000 stochastic samples of different random networks for each edge density. For parameters, β = 0.5, g 0 = 2, d 0 = −1, ⟨ g ⟩ = −0.25. Spatial drug heterogeneity effects and optimal strategy under a given multi-microhabitat structure Different spatial drug heterogeneities, like different multi-microhabitat structures, induce distinct population responses. In our study, to allow meaningful comparison across scenarios, we control the spatially averaged drug dose - or equivalently, the spatially averaged growth rate ⟨ g ⟩ - in our minimal model. We consider two cases depending on whether spatial drug heterogeneity is controllable: 1. Fully controllable heterogeneity: Clinically, if we are fortunate enough to precisely control drug concentration at each microhabitat (or node), population decline can be achieved by selecting the spatial drug assignment that minimizes the largest eigenvalue - we simply find when min λ max ( A ) < 1. This gives us the optimal spatial drug configuration among all possible combinations. 2. Uncontrollable or unknown heterogeneity: If spatial drug assignment cannot be precisely controlled or is unknown, we must ensure population decline under the worst-case spatial drug configuration. This requires max λ max ( A ) < 1, ensuring a robust population decline that is independent of the specific drug distribution. This naturally leads to a constrained optimization problem, as illustrated in recent work on robust decline under spatial drug heterogeneity in continuous 1D space [ 44 ]. In both cases, we are solving a constrained nonlinear optimization problem involving the largest eigenvalue. However, the search space is large, since each g i varies continuously and the configuration space grows combinatorially with the number of microhabitats N . Remarkably, the solution to this class of problems can be simplified. For the minimization case (min λ max ( A )), assuming ⟨ g ⟩ < 0, the optimal strategy is alway an interior point with every 0 < δ i < c , or d 0 < g i < g 0 . This optimal strategy depends on K and is usually hard to derive analytically. However, a uniform drug assignment: g i = ⟨ g ⟩ for all microhabitats, as the “even-spread strategy”, is good enough and can always induce population decline under ⟨ g ⟩ < 0. For the maximization case (max λ max ( A )), the optimum typically lies on or near the active constraints (see SI for a detailed proof). This means that, for a fixed ⟨ g ⟩, we should assign as many microhabitats as possible to the maximum death rate d 0 or the maximum growth rate g 0 , leaving at most one remainder g r ∈ [ d 0 , g 0 ]. The optimal configuration then takes the form { g i } optimal = { g r ; g 0 …, …; d 0 ,}, which we refer to as the “one-remainder strategy”. Hence, this continuous eigenvalue optimization problem reduces to a discrete combinatorial optimization. Interestingly, our kernel transformation is able to further simplify the problem. By expressing growth in terms of relative rates δ i , the optimal growth configuration becomes { δ i } optimal = { δ r ; c ,…; 0,…}, where δ r ∈[0, c ], and is the maximum possible relative growth rate. For nodes assigned g i = d 0 , δ i = 0, which effectively removes all edges from microhabitat i to its neighbors under our transformation (as the structure becomes effectively fully connected). If N − S microhabitats are suppressed ( δ i = 0), the remaining effective graph forms an S -clique, and the matrix A reduces to a smaller principal submatrix A S of size S . This matrix can be decomposed as: Here, we assume the microhabitat with the remainder g r is labeled , and k 1 = ( K 11 , K 21 , …, K S 1 ) T is the first column of K [ S ]. K 11 reflects the inverse dynamic-related centrality of the node with the remainder. Thus, solving the optimization problem becomes equivalent to selecting an S -clique from the full graph to minimize or maximize λ max ( A S ). While discrete combinatorics optimization remains computationally challenging, we propose a centrality-based heuristic that leverages the decoupling of structure and growth. Under very high drug doses - i.e., when , or when g 0 ≥ ( N − 1) | d 0 | - the one-remainder strategy becomes equivalent to selecting a single growing microhabitat. Then we have λ max ( A ) = δ r K ii . To guarantee decline under unknown heterogeneity, we require , implying the structure must contain at least one relatively isolated microhabitat. In the general case with multiple positive growth rates { δ r ; c ,…}, we aim to identify an interconnected but relatively isolated S -clique to ensure robust decline. This intuition aligns with Figure 3B , where less connected graphs tend to suppress growth less effectively, thus serving as the “worst” case. Note that selecting the most isolated S-clique is not simply equivalent to minimizing the row sum of A S , which is only an upper bound: . Optimal drug assignment with incomplete controllable microhabitats Under complete control, the good enough strategy to minimize λ max ( A ) is to assign a uniform drug distribution: g i = ⟨ g ⟩ across all microhabitats. However, in practice, we may only be able to control a subset of microhabitats and apply high enough drug concentrations to push them to the maximum death rate d 0 , while the remaining microhabitats stay approximately drug-free at growth rate g 0 . In this scenario, the problem again becomes a discrete combinatorial optimization - selecting the optimal subcore. When the number of controllable microhabitats is large, we should avoid targeting the well-connected core. Instead, we preserve this core (i.e., assign it positive growth rates), since it can effectively “sense” the largest number of the surrounding “sink” microhabitats with maximum death rate d 0 . Interestingly, when the number of controllable microhabitats is very limited-e.g., we can apply high drug concentration to only one node (setting δ i = 0)-we find that the optimal strategy is to suppress the microhabitat with the smallest K ii , i.e., the node with the highest dynamic-related centrality (which often corresponds to highest degree). With more available microhabitats to control, we continue suppressing other nodes with high centralities in descending order. Table 1 summarizes optimized strategies across different graph families based on this principle. For more complex structures, this centrality-based rule often provides a greedy but effective approximation of the optimal solution. However, to find the true optimal drug configuration, one still needs to solve the full discrete combinatorial problem of minimizing λ ( A S ) with partial control. View this table: View inline View popup Download powerpoint Table 1. Optimized clearance strategies with limited controllable microhabitats, across different graph families based on the regularized Laplacian kernel matrix K This drug assignment strategy under limited controllability closely mirrors strategies used in epidemic outbreak suppression in susceptible-infectious-susceptible (SIS) network dynamics [ 70 – 73 ], where the optimal curing rate of each node is made proportional to its degree. Thus, in our minimal growth-migration model: a. when many microhabitats are controllable, we protect the well-connected core; b. when only a few are controllable, we target the core to suppress it most effectively. Robust population decline with unknown spatial drug heterogeneity In most real-world scenarios, especially within the human body, estimating or controlling drug concentrations at each microhabitat is infeasible. Spatial drug heterogeneity is typically unknown ( Figure 4A ). Assuming we can estimate the migration rate β and control the total drug dose (or the spatially averaged growth rate ⟨ g ⟩), the problem becomes identifying robust clearance strategies that ensure population decline in the worst case. That is, we need: Download figure Open in new tab Figure 4. Robust population decline under unknown spatial drug heterogeneity. A. Illustration of the concept of robust population decline. When spatial drug heterogeneity is unknown, the clearance strategy evaluates whether population decline always occurs under a given migration rate and total drug dose, effectively eliminating the variability induced by different spatial drug configurations. B . Example phase diagram on a cycle graph, showing the relationship between migration rate β and spatially averaged growth rate ⟨ g ⟩. The bottom-right corner presents the mathematical form of the sufficient condition, which depends on β , maximum death rate d 0 , minimum degree k min , maximum Laplacian eigenvalue ω max , and maximum relative growth rate δ max . Above the phase diagram, four final population densities after a finite time are shown for different spatial growth distributions selected from robust decline phase. The dashed line indicates the initial uniform population density. On the right, the corresponding spatial growth rate distributions are illustrated. C . Phase diagrams for two additional structures — star and complete graphs — under relatively high maximum growth rate g 0 = 2( N −1) and low maximum death rate d 0 = −1. For the complete graph, the sufficient condition becomes exact, perfectly matching the numerical boundary. This is again a constrained optimization problem, which, as discussed earlier, reduces to selecting highly interconnected subcores (see SI for complete proof). To illustrate this framework across multi-microhabitat structures, we present migration-growth phase diagrams ( β vs. ⟨ g ⟩) that identify parameter regimes yielding robust population decline under ⟨ g ⟩ < 0. Besides robust population decline, population responses may still vary depending on heterogeneity, resulting in a “mixed phase” region where both growth and decline are possible. Note that for any ⟨ g ⟩ < 0, the “even-spread” strategy always achieves decline, which means there’s no region where robust population growth is possible. For simplicity, we focus on the case where g 0 ≥( N −1) | d 0 |, or equivalently c ≥ N . In this regime, the criterion: becomes an exact condition for robust population decline. Here, δ r is the maximum possible relative growth rate, and max i K ii is the inverse dynamic-related centrality of the most isolated microhabitat. This condition tells us that as long as max i K ii < 1 /δ r , we are guaranteed robust population decline. Since K ii ∝ β −1 and δ r ∝ ⟨ g ⟩, this condition reveals a tradeoff: higher migration suppresses growth, while higher average growth promotes persistence. To better visualize this interplay, we apply an interpolation inequality to max i K ii and rearrange terms (see SI), yielding a simplified analytical boundary for robust decline based only on structural extremes: Here, k min is the minimum degree and ω max the largest Laplacian eigenvalue of the structure. We explicitly write δ max = δ r to emphasize its role as the extreme growth parameter. This bound clearly captures the fundamental tradeoff: to suppress populations with high intrinsic growth ( δ r ), we need stronger inter-microhabitat migration ( β ). For a complete graph, where all nonzero Laplacian eigenvalues are equal, this inequality becomes exact (see SI). In Figure 4B , we illustrate the phase diagram for a 12-node cycle graph. The blue region (top left) corresponds to robust decline (max λ max ( A ) < 1), while the white region represents the mixed phase. Our analytical bound provides a higher estimate for the true numerical boundary. In cancer research, the star graph is especially relevant for modeling metastasis from a primary site to distant tissues [ 49 , 51 ]. Similarly, in bacterial communities, fully connected networks may arise [ 74 , 75 ]. Figure 4C presents phase diagrams for both star and complete graphs. In the complete graph (bottom panel), the sufficient condition precisely matches the numerical boundary. For these two structures, we can derive closed-form expressions for max i K ii . For star graph, the largest K is at a leaf, and , where . For complete graph, all nodes are identical, and . Thus, for graphs like these, the boundary line can be exactly determined by solving δ r · max i K ii = 1. For more complex graphs where max i K ii is not analytically tractable (e.g., even cycles), our extreme-bound approximation from equation (8) still provides a fast and reliable method to identify regions of robust population decline. Discussion Complex biological systems, such as cancer and bacterial populations, are deeply intertwined with network structure, growth dynamics, and spatial interactions. In this study, we introduced a kernel transformation that enabled the derivation of a new theoretical criterion predicting population-level responses across arbitrary multi-microhabitat structures under spatially heterogeneous drug environments. This transformation naturally links population dynamics to a centrality measure known as forest closeness centrality, which we refer to as dynamic-related centrality in our framework. Our findings demonstrate that increasing connectivity between microhabitats consistently promotes population decline, effectively mimicking the impact of increasing the migration rate. Furthermore, we show that the optimization of spatial drug heterogeneity for inducing population decline can be reduced to a subgraph selection problem on the transformed fully connected graph. In addition, we identify a parameter regime in the growth-migration phase diagram that supports robust population decline — a guaranteed clearance across all spatial drug configurations if the actual spatial drug heterogeneity is unknown. We derive a sufficient theoretical condition characterizing this phase, capturing the interplay between spatially averaged growth and migration rates. Notably, this condition becomes exact when the underlying microhabitat structure is a complete graph. In [ 44 ], perturbation approximation was successfully applied to derive an analytical expression for the largest eigenvalue λ max ( G − βL ), where is the Laplacian operator with two absorbing boundaries on a one-dimensional continuous space. The approximated form revealed two optimized strategies for spatial drug arrangement by examining the unperturbed largest eigenvector of L : the “importance” of spatial positions is ranked by the squared components of this eigenvector. And it’s found that, assigning drug-free growth rates g 0 to central positions tends to mitigate population decline, while assigning g 0 near absorbing boundaries tends to accelerate it. However, when the spatial structure becomes discrete, and the driving force of decline arises not from absorbing boundaries, but from “sink” microhabitats with negative growth rates, perturbation theory fails to accurately capture the system’s behavior—it reflects only the effect of spatially averaged growth. In contrast, our kernel transformation introduces a new matrix A = D 1 / 2 KD 1 / 2 , where K = ( I + γL ) −1 is the regularized Laplacian kernel and D is the diagonal matrix of growth-related weights. Even first-order perturbation theory fails on K , since K and L share the same eigenbasis, and thus the largest eigenvector still assigns equal importance to all microhabitats. Despite this, the diagonal elements K ii , which correspond to forest closeness centrality — or what we refer to as “dynamic-related centrality” — provide a meaningful measure for ranking the relative importance of microhabitats. Although perturbation theory is insufficient in this setting, the kernel matrix K still allows us to heuristically infer microhabitat importance and guide spatial drug dosing strategies. In particular, it enables a greedy subcore selection strategy on the transformed fully connected graph, as discussed in the main text. This insight allows us to extend the one-dimensional intuition: assigning g 0 to central or boundary positions depends on the controllability context. In our discrete framework, the “center” corresponds to a well-interconnected core, while the “edges” represent relatively isolated subgraphs. Thus, through kernel transformation, we heuristically recover spatial clearance strategies based on dynamic-related centralities, even in cases where classical perturbation methods break down. Across-node properties such as centrality measures play a crucial role in graph theory and network analysis. For instance, subgraph centrality quantifies node importance based on participation in subgraphs [ 76 ], with a bias toward smaller subgraphs. Interestingly, our dynamic-related centrality C exhibits a similar relationship with traditional closeness centrality, correlating with node degree, and other centrality measures such as Katz or eigenvalue centrality [ 77 , 78 ]. Recent work on dynamic centrality has aimed to integrate both topological and dynamical properties of nodes to more accurately identify influential spreaders in static and temporal networks [ 79 – 81 ]. These studies emphasize the importance of incorporating temporal aspects into centrality metrics, as static measures fail to capture the evolving nature of real-world networks, crucial for predicting node importance in dynamic systems [ 80 , 82 ]. Empirical studies demonstrate that these methods outperform conventional metrics across various scenarios, including information and disease spreading [ 81 , 82 ]. However, previous formulations of dynamic centrality primarily focus on the addition or deletion of nodes within a network, meaning they remain fundamentally structural measures that describe the dynamics of the network itself. In contrast, our dynamic-related centrality focuses on predicting the results of dynamics occurring on the network, and it naturally emerges from the dynamics itself by doing the kernel transformation. This distinction makes our proposed measure novel, and we term it “dynamic-related centrality” to distinguish it from “dynamic centrality” for the dynamics of the network. We anticipate that this new centrality will expand existing centrality classifications and provide direct insights into dynamic properties without requiring complex simulations. Interestingly, our matrix K , which defines the dynamic-related centrality C , shares the same mathematical form as the classical regularized Laplacian kernel (RLK) [ 83 ] — a well-known graph similarity measure based on the Laplacian matrix L . The RLK is defined as K = ( I + γL ) −1 , where I is the identity matrix, L is the Laplacian, and γ > 0 is a tunable regularization parameter that controls the global influence of the Laplacian. Despite this shared mathematical form, we emphasize important conceptual differences. First, in our growth-migration framework, the parameter arises naturally from the underlying population dynamics, rather than being introduced as a free parameter. Second, while the RLK is typically used to analyze off-diagonal elements to quantify node-to-node similarity, our analysis centers on the diagonal entries K ii , which correspond to forest closeness centrality and reflect the self-centrality of each node in the dynamic context (While λ max ( A ) can be interpreted through a walk-based perspective—capturing contributions from paths traversing multiple nodes thus considering interactions between different nodes — this interpretation is beyond the scope of the present work and is not central to our theoretical framework). Establishing these connections between our dynamic-related centrality and classical graph-theoretic measures not only enhances the interpretability of our framework but also expands its potential applications in network-based modeling of biological systems. In the main text, we did not emphasize the upper bound , which not only serves as a useful heuristic for guiding optimized drug allocation strategies, but also reveals deeper structural insights. Specifically, we find that this maximum row sum is tightly linked to the presence of local symmetries in the system. When the matrix A , or a principal submatrix A [ S ], is circulant, the inequality becomes an equality: . Although spatial drug heterogeneity generally breaks most global symmetries in the underlying structure, certain local symmetries may still be preserved. In such cases, this equality offers a convenient and accurate method for predicting population response by directly evaluating the maximum row sum of A . Detailed derivations, along with illustrative examples, are provided in the Supplementary Information. Understanding population responses to drug treatment, such as antibiotics or chemotherapy, is critical for assessing therapeutic success. Recent studies in both bacterial and cancer systems have identified an interesting phenomenon known as “bistability”, where small differences in key parameters—such as drug concentration or growth rate—can lead to drastically different outcomes [ 84 – 88 ]. This concept helps explain variability in clinical outcomes. Here, we propose an alternative explanation for such bi-outcome scenarios: spatial drug heterogeneity and differences in microhabitat connectivity may also contribute to distinct clearance outcomes, independent of bistability. This insight could inform the development of new treatment strategies tailored to spatial drug distribution patterns. Experimental studies have demonstrated the impact of spatial structure on growth [ 44 , 89 ], competition, and evolution [ 90 – 93 ]. Our theoretical findings may aid in the interpretation of experimental results and guide the design of future studies. Combining theoretical and experimental approaches can enhance our understanding and contribute to the development of spatially explicit drug dosing strategies. This study explored population responses across various microhabitat structures and spatial drug heterogeneities, revealing insights beyond previous research that primarily focused on homogeneous structures and drug distributions. However, several extensions remain for future research. First, we considered only non-frequency-dependent drug-induced growth in a single cell type. Although experimental or clinical data on drug-mediated growth interactions remain scarce, some recent studies have begun to address this gap [ 94 , 95 ]. Future work could incorporate frequency-dependent growth effects and multi-species interactions using evolutionary game theory [ 96 – 100 ] to uncover additional theoretical results and experimental comparisons. Second, we did not account for temporal drug decay or fluctuations in our models. Future work could integrate spatial-temporal drug variations based on patient-specific pharmacokinetics/pharmacodynamics (PK/PD) [ 25 ] to inform personalized treatment strategies that ensure bacterial or cancer clearance. Acknowledgments This study was supported by NIH R35GM124875 (K. W.). We’d like to sincerely thank suggestions from Dr. David Lubensky and Dr. Michal Zochowoski. Footnotes Corrected calculated errors in major findings ↵ 1 , with V and {ω k } representing the eigenvector matrix and eigenvalues of the Laplacian matrix L , respectively. It contains only structure information, with maximum drug-induced death rate | d 0 | as the baseline growth. References 1. ↵ K. Wood , “Microbial ecology: complex bacterial communities reduce selection for antibiotic resistance ,” Current biology , vol. 29 , no. 21 , pp. R1143 – R1145 , 2019 . OpenUrl CrossRef PubMed 2. M. Denk-Lobnig and K. B. Wood , “Antibiotic resistance in bacterial communities ,” Current Opinion in Microbiology , vol. 74 , p. 102306 , 2023 . OpenUrl CrossRef PubMed 3. ↵ G. Aguadé-GorgoriÓ , A. R. Anderson , and R. Solé , “Modeling tumors as complex ecosystems ,” Iscience , vol. 27 , no. 9 , 2024 . 4. ↵ E. Hansen , J. Karslake , R. J. Woods , A. F. Read , and K. B. Wood , “Antibiotics can be used to contain drug-resistant bacteria by maintaining sufficiently large sensitive populations ,” PLoS biology , vol. 18 , no. 5 , p. e3000713 , 2020 . OpenUrl CrossRef PubMed 5. ↵ J. A. Gallaher , P. M. Enriquez-Navas , K. A. Luddy , R. A. Gatenby , and A. R. Anderson , “Spatial heterogeneity and evolutionary dynamics modulate time to recurrence in continuous and adaptive cancer therapies ,” Cancer research , vol. 78 , no. 8 , pp. 2127 – 2139 , 2018 . OpenUrl Abstract / FREE Full Text 6. J. West , F. Adler , J. Gallaher , M. Strobl , R. Brady-Nicholls , J. Brown , M. Roberson-Tessi , E. Kim , R. Noble , Y. Viossat , et al. , “A survey of open questions in adaptive therapy: Bridging mathematics and clinical translation ,” Elife , vol. 12 , p. e84263 , 2023 . OpenUrl CrossRef PubMed 7. M. A. Strobl , A. L. Martin , J. West , J. Gallaher , M. Robertson-Tessi , R. Gatenby , R. Wenham , P. K. Maini , M. Damaghi , and A. R. Anderson , “To modulate or to skip: De-escalating parp inhibitor maintenance therapy in ovarian cancer using adaptive therapy ,” Cell Systems , vol. 15 , no. 6 , pp. 510 – 525 , 2024 . OpenUrl CrossRef PubMed 8. ↵ K. Gallagher , M. A. Strobl , D. S. Park , F. C. Spoendlin , R. A. Gatenby , P. K. Maini , and A. R. Anderson , “Mathematical model-driven deep learning enables personalized adaptive therapy ,” Cancer Research , vol. 84 , no. 11 , pp. 1929 – 1941 , 2024 . OpenUrl CrossRef PubMed 9. ↵ P. P. Chakraborty , L. R. Nemzer , and R. Kassen , “Experimental evidence that metapopulation structure can accelerate adaptive evolution ,” BioRxiv , pp. 2021 – 07 , 2021 . 10. ↵ A. Wu , K. Loutherback , G. Lambert , L. Estévez-SalmerÓn , T. D. Tlsty , R. H. Austin , and J. C. Sturm , “Cell motility and drug gradients in the emergence of resistance to chemotherapy ,” Proceedings of the National Academy of Sciences , vol. 110 , no. 40 , pp. 16103 – 16108 , 2013 . OpenUrl Abstract / FREE Full Text 11. ↵ Q. Zhang , G. Lambert , D. Liao , H. Kim , K. Robin , C.-k. Tung , N. Pourmand , and R. H. Austin , “Acceleration of emergence of bacterial antibiotic resistance in connected microenvironments ,” Science , vol. 333 , no. 6050 , pp. 1764 – 1767 , 2011 . OpenUrl Abstract / FREE Full Text 12. ↵ M. Baym , T. D. Lieberman , E. D. Kelsic , R. Chait , R. Gross , I. Yelin , and R. Kishony , “Spatiotemporal microbial evolution on antibiotic landscapes ,” Science , vol. 353 , no. 6304 , pp. 1147 – 1151 , 2016 . OpenUrl Abstract / FREE Full Text 13. K. Nagy , B. Dukic , O. Hodula , Á. Ábrahám , E. Csákvári , L. Dér , M. T. Wetherington , J. Noorlag , J. E. Keymer , and P. Galajda , “Emergence of resistant escherichia coli mutants in microfluidic on-chip antibiotic gradients ,” Frontiers in microbiology , vol. 13 , p. 820738 , 2022 . OpenUrl CrossRef PubMed 14. G. G. Perron , A. Gonzalez , and A. Buckling , “Source–sink dynamics shape the evolution of antibiotic resistance and its pleiotropic fitness cost ,” Proceedings of the Royal Society B: Biological Sciences , vol. 274 , no. 1623 , pp. 2351 – 2356 , 2007 . OpenUrl CrossRef PubMed Web of Science 15. ↵ R. M. Wheatley , J. D. Caballero , T. E. van der Schalk , F. H. De Winter , L. P. Shaw , N. Kapel , C. Recanatini , L. Timbermont , J. Kluytmans , M. Esser , et al. , “Gut to lung translocation and antibiotic mediated selection shape the dynamics of pseudomonas aeruginosa in an icu patient ,” Nature communications , vol. 13 , no. 1 , p. 6523 , 2022 . OpenUrl CrossRef PubMed 16. ↵ R. Hermsen and T. Hwa , “Sources and sinks: a stochastic model of evolution in heterogeneous environments ,” Physical review letters , vol. 105 , no. 24 , p. 248104 , 2010 . OpenUrl CrossRef PubMed 17. ↵ R. Hermsen , J. B. Deris , and T. Hwa , “On the rapidity of antibiotic resistance evolution facilitated by a concentration gradient ,” Proceedings of the National Academy of Sciences , vol. 109 , no. 27 , pp. 10775 – 10780 , 2012 . OpenUrl Abstract / FREE Full Text 18. D. Labavić , C. Loverdo , and A.-F. Bitbol , “Hydrodynamic flow and concentration gradients in the gut enhance neutral bacterial diversity ,” Proceedings of the National Academy of Sciences , vol. 119 , no. 1 , p. e2108671119 , 2022 . OpenUrl CrossRef PubMed 19. F. Baquero , M.-C. Negri , M.-I. Morosini , and J. Blázquez , “Antibiotic-selective environments ,” Clinical infectious diseases , vol. 27 , no. Supplement 1, pp. S5 – S11 , 1998 . OpenUrl CrossRef PubMed 20. T. B. Kepler and A. S. Perelson , “Drug concentration heterogeneity facilitates the evolution of drug resistance ,” Proceedings of the National Academy of Sciences , vol. 95 , no. 20 , pp. 11514 – 11519 , 1998 . OpenUrl Abstract / FREE Full Text 21. P. Greulich , B. Waclaw , and R. J. Allen , “Mutational pathway determines whether drug gradients accelerate evolution of drug-resistant cells ,” Physical review letters , vol. 109 , no. 8 , p. 088101 , 2012 . OpenUrl CrossRef PubMed 22. L. Marrec and A.-F. Bitbol , “Adapt or perish: Evolutionary rescue in a gradually deteriorating environment ,” Genetics , vol. 216 , no. 2 , pp. 573 – 583 , 2020 . OpenUrl Abstract / FREE Full Text 23. H. Steel and A. Papachristodoulou , “The effect of spatiotemporal antibiotic inhomogeneities on the evolution of resistance ,” Journal of theoretical biology , vol. 486 , p. 110077 , 2020 . OpenUrl CrossRef PubMed 24. A. I. Klimenko , Y. G. Matushkin , N. A. Kolchanov , and S. A. Lashin , “Spatial heterogeneity promotes antagonistic evolutionary scenarios in microbial community explained by ecological stratification: a simulation study ,” Ecological modelling , vol. 399 , pp. 66 – 76 , 2019 . OpenUrl CrossRef 25. ↵ H. Xie , Y. Jiao , Q. Fan , M. Hai , J. Yang , Z. Hu , Y. Yang , J. Shuai , G. Chen , R. Liu , et al. , “Modeling three-dimensional invasive solid tumor growth in heterogeneous microenvironment under chemotherapy ,” PloS one , vol. 13 , no. 10 , p. e0206292 , 2018 . OpenUrl CrossRef PubMed 26. J. Deng , L. Zhou , R. A. Sanford , L. A. Shechtman , Y. Dong , R. E. Alcalde , M. Sivaguru , G. A. Fried , C. J. Werth , and B. W. Fouke , “Adaptive evolution of escherichia coli to ciprofloxacin in controlled stress environments: contrasting patterns of resistance in spatially varying versus uniformly mixed concentration conditions ,” Environmental Science & Technology , vol. 53 , no. 14 , pp. 7996 – 8005 , 2019 . OpenUrl CrossRef PubMed 27. A. M. Leale and R. Kassen , “The emergence, maintenance, and demise of diversity in a spatially variable antibiotic regime ,” Evolution Letters , vol. 2 , no. 2 , pp. 134 – 143 , 2018 . OpenUrl CrossRef PubMed 28. R. C. MacLean and T. Vogwill , “Limits to compensatory adaptation and the persistence of antibiotic resistance in pathogenic bacteria ,” Evolution, medicine, and public Health , vol. 2015 , no. 1 , pp. 4 – 12 , 2015 . OpenUrl CrossRef PubMed 29. S. Xu , Evolution of Microbial Interactions . PhD thesis, University of Miami , 2020 . 30. B. Zhang , L. Zhai , J. Bintz , S. M. Lenhart , W. Valega-Mackenzie , and J. D. Van Dyken , “The optimal controlling strategy on a dispersing population in a two-patch system: Experimental and theoretical perspectives ,” Journal of theoretical biology , vol. 528 , p. 110835 , 2021 . OpenUrl CrossRef PubMed 31. A. Abbara and A.-F. Bitbol , “Frequent asymmetric migrations suppress natural selection in spatially structured populations ,” bioRxiv , pp. 2023 – 06 , 2023 . 32. P. Amarasekare , “Competitive coexistence in spatially structured environments: a synthesis ,” Ecology letters , vol. 6 , no. 12 , pp. 1109 – 1122 , 2003 . OpenUrl CrossRef Web of Science 33. Y. Liang , S. McDonnell , and M. Clynes , “Examining the relationship between cancer invasion/metastasis and drug resistance ,” Current cancer drug targets , vol. 2 , no. 3 , pp. 257 – 277 , 2002 . OpenUrl CrossRef PubMed 34. J. Hansen , Heterogeneity and the Evolution of Drug Resistance . PhD thesis, University of Michigan , 2023 . 35. F. Fu , M. A. Nowak , and S. Bonhoeffer , “Spatial heterogeneity in drug concentrations can facilitate the emergence of resistance to cancer therapy ,” PLoS computational biology , vol. 11 , no. 3 , p. e1004142 , 2015 . OpenUrl CrossRef 36. ↵ V. Piskovsky and N. M. Oliveira , “Bacterial motility can govern the dynamics of antibiotic resistance evolution ,” Nature Communications , vol. 14 , no. 1 , p. 5584 , 2023 . OpenUrl CrossRef PubMed 37. L. Marrec , I. Lamberti , and A.-F. Bitbol , “Toward a universal model for spatially structured populations ,” Physical review letters , vol. 127 , no. 21 , p. 218102 , 2021 . OpenUrl CrossRef PubMed 38. O. M. Ghosh and B. H. Good , “Emergent evolutionary forces in spatial models of luminal growth and their application to the human gut microbiota ,” Proceedings of the National Academy of Sciences , vol. 119 , no. 28 , p. e2114931119 , 2022 . OpenUrl CrossRef PubMed 39. ↵ M. G. De Jong and K. B. Wood , “Tuning spatial profiles of selection pressure to modulate the evolution of drug resistance ,” Physical review letters , vol. 120 , no. 23 , p. 238102 , 2018 . OpenUrl CrossRef PubMed 40. ↵ M. De Jong , Emergence of Antibiotic Resistance in Microbial Populations with Spatial Heterogeneity . PhD thesis, University of Michigan , 2021 . 41. ↵ J. T. Lamont , “Patient education: Helicobacter pylori infection and treatment (beyond the basics),” 2022 . 42. ↵ R. Dent , M. Trudeau , K. I. Pritchard , W. M. Hanna , H. K. Kahn , C. A. Sawka , L. A. Lickley , E. Rawlinson , P. Sun , and S. A. Narod , “Triple-negative breast cancer: clinical features and patterns of recurrence ,” Clinical cancer research , vol. 13 , no. 15 , pp. 4429 – 4434 , 2007 . OpenUrl Abstract / FREE Full Text 43. ↵ A. I. Riggio , K. E. Varley , and A. L. Welm , “The lingering mysteries of metastatic recurrence in breast cancer ,” British journal of cancer , vol. 124 , no. 1 , pp. 13 – 26 , 2021 . OpenUrl CrossRef PubMed 44. ↵ Z. Hu , Y. Wu , T. Freire , E. Gjini , and K. Wood , “Linking spatial drug heterogeneity to microbial growth dynamics in theory and experiment ,” bioRxiv , pp. 2024 – 11 , 2024 . 45. ↵ T. F. A. Freire , Z. Hu , K. B. Wood , and E. Gjini , “Modeling spatial evolution of multi-drug resistance under drug environmental gradients ,” PLOS Computational Biology , vol. 20 , no. 5 , p. e1012098 , 2024 . OpenUrl CrossRef PubMed 46. ↵ L. Chen , N. Blumm , N. Christakis , A. Barabasi , and T. S. Deisboeck , “Cancer metastasis networks and the prediction of progression patterns ,” British journal of cancer , vol. 101 , no. 5 , pp. 749 – 758 , 2009 . OpenUrl CrossRef PubMed Web of Science 47. J. Budczies , M. von Winterfeld , F. Klauschen , M. Bockmayr , J. K. Lennerz , C. Denkert , T. Wolf , A. Warth , M. Dietel , I. Anagnostopoulos , et al. , “The landscape of metastatic progression patterns across major human cancers ,” Oncotarget , vol. 6 , no. 1 , p. 570 , 2014 . OpenUrl 48. L. d. F. Costa , O. N. Oliveira Jr , G. Travieso , F. A. Rodrigues , P. R. Villas Boas , L. Antiqueira , M. P. Viana , and L. E. Correa Rocha , “Analyzing and modeling real-world phenomena with complex networks: a survey of applications ,” Advances in Physics , vol. 60 , no. 3 , pp. 329 – 412 , 2011 . OpenUrl CrossRef 49. ↵ M. Riihimäki , A. Hemminki , M. Fallah , H. Thomsen , K. Sundquist , J. Sundquist , and K. Hemminki , “Metastatic sites and survival in lung cancer ,” Lung cancer , vol. 86 , no. 1 , pp. 78 – 84 , 2014 . OpenUrl CrossRef PubMed 50. M. Riihimäki , H. Thomsen , K. Sundquist , J. Sundquist , and K. Hemminki , “Clinical landscape of cancer metastases ,” Cancer medicine , vol. 7 , no. 11 , pp. 5534 – 5542 , 2018 . OpenUrl CrossRef PubMed 51. ↵ X. Jin , Z. Demere , K. Nair , A. Ali , G. B. Ferraro , T. Natoli , A. Deik , L. Petronio , A. A. Tang , C. Zhu , et al. , “A metastasis map of human cancer cell lines ,” Nature , vol. 588 , no. 7837 , pp. 331 – 336 , 2020 . OpenUrl CrossRef PubMed 52. ↵ N. L. Trevaskis , L. M. Kaminskas , and C. J. Porter , “From sewer to saviour—targeting the lymphatic system to promote drug exposure and activity ,” Nature reviews Drug discovery , vol. 14 , no. 11 , pp. 781 – 803 , 2015 . OpenUrl CrossRef PubMed 53. T. Kodama , D. Matsuki , A. Tada , K. Takeda , and S. Mori , “New concept for the prevention and treatment of metastatic lymph nodes using chemotherapy administered via the lymphatic network ,” Scientific reports , vol. 6 , no. 1 , p. 32506 , 2016 . OpenUrl CrossRef PubMed 54. ↵ H. Ji , C. Hu , X. Yang , Y. Liu , G. Ji , S. Ge , X. Wang , and M. Wang , “Lymph node metastasis in cancer progression: molecular mechanisms, clinical significance and therapeutic interventions ,” Signal Transduction and Targeted Therapy , vol. 8 , no. 1 , p. 367 , 2023 . OpenUrl CrossRef 55. ↵ M. Ziaka and A. Exadaktylos , “Exploring the lung-gut direction of the gut-lung axis in patients with ards ,” Critical Care , vol. 28 , no. 1 , p. 179 , 2024 . OpenUrl CrossRef PubMed 56. ↵ M. W. Eladham , B. Selvakumar , N. S. Sharif-Askari , F. S. Sharif-Askari , S. M. Ibrahim , and R. Halwani , “Unraveling the gut-lung axis: Exploring complex mechanisms in disease interplay ,” Heliyon , 2024 . 57. ↵ J. Coates , B. R. Park , D. Le , E. Şimşek , W. Chaudhry , and M. Kim , “Antibiotic-induced population fluctuations and stochastic clearance of bacteria ,” Elife , vol. 7 , p. e32976 , 2018 . OpenUrl CrossRef PubMed 58. ↵ N. G. Van Kampen , Stochastic processes in physics and chemistry , vol. 1 . Elsevier , 1992 . 59. ↵ A. J. Lee , S. Wang , H. R. Meredith , B. Zhuang , Z. Dai , and L. You , “Robust, linear correlations between growth rates and β-lactam–mediated lysis rates ,” Proceedings of the National Academy of Sciences , vol. 115 , no. 16 , pp. 4069 – 4074 , 2018 . OpenUrl Abstract / FREE Full Text 60. S. Pearl Mizrahi , A. Goyal , and J. Gore , “Community interactions drive the evolution of antibiotic tolerance in bacteria ,” Proceedings of the National Academy of Sciences , vol. 120 , no. 3 , p. e2209043119 , 2023 . OpenUrl CrossRef PubMed 61. ↵ A. Bren , D. S. Glass , Y. K. Kohanim , A. Mayo , and U. Alon , “Tradeoffs in bacterial physiology determine the efficiency of antibiotic killing ,” Proceedings of the National Academy of Sciences , vol. 120 , no. 51 , p. e2312651120 , 2023 . OpenUrl CrossRef PubMed 62. ↵ L. Geyrhofer , P. Ruelens , A. D. Farr , D. Pesce , J. A. G. de Visser , and N. Brenner , “Minimal surviving inoculum in collective antibiotic resistance ,” Mbio , vol. 14 , no. 2 , pp. e02456 – 22 , 2023 . OpenUrl PubMed 63. ↵ K. Quinn , E. Terzi , and M. Crovella , “ Characterizing covid waves via spatio-temporal decomposition ,” in Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining , pp. 3783 – 3791 , 2022 . 64. ↵ F. Aziz , L. T. Slater , L. Bravo-Merodio , A. Acharjee , and G. V. Gkoutos , “Link prediction in complex network using information flow ,” Scientific Reports , vol. 13 , no. 1 , p. 14660 , 2023 . OpenUrl CrossRef PubMed 65. ↵ L. Huang , L. Liao , and C. H. Wu , “Protein-protein interaction prediction based on multiple kernels and partial network with linear programming ,” BMC Systems Biology , vol. 10 , pp. 151 – 164 , 2016 . OpenUrl 66. ↵ J. Fan , A. Cannistra , I. Fried , T. Lim , T. Schaffner , M. Crovella , B. Hescott , and M. D. Leiserson , “Functional protein representations from biological networks enable diverse cross-species inference ,” Nucleic acids research , vol. 47 , no. 9 , pp. e51 – e51 , 2019 . OpenUrl CrossRef PubMed 67. ↵ S. Ray , S. Lall , and S. Bandyopadhyay , “A deep integrated framework for predicting sars-cov2–human protein-protein interaction ,” IEEE Transactions on Emerging Topics in Computational Intelligence , vol. 6 , no. 6 , pp. 1463 – 1472 , 2022 . OpenUrl 68. ↵ Y. Jin , Q. Bao , and Z. Zhang , “ Forest distance closeness centrality in disconnected graphs ,” in 2019 IEEE international conference on data mining (icdm) , pp. 339 – 348 , IEEE , 2019 . 69. ↵ A. van der Grinten , E. Angriman , M. Predari , and H. Meyerhenke , “ New approximation algorithms for forest closeness centrality–for individual vertices and vertex groups ,” in Proceedings of the 2021 SIAM International Conference on Data Mining (SDM) , pp. 136 – 144 , SIAM , 2021 . 70. ↵ S. C. Ferreira , C. Castellano , and R. Pastor-Satorras , “Epidemic thresholds of the susceptible-infected-susceptible model on networks: A comparison of numerical and theoretical results ,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics , vol. 86 , no. 4 , p. 041125 , 2012 . OpenUrl 71. C. Li , R. van de Bovenkamp , and P. Van Mieghem , “Susceptible-infected-susceptible model: A comparison of n-intertwined and heterogeneous mean-field approximations ,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics , vol. 86 , no. 2 , p. 026116 , 2012 . OpenUrl 72. E. Cator and P. Van Mieghem , “Susceptible-infected-susceptible epidemics on the complete graph and the star graph: Exact analysis ,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics , vol. 87 , no. 1 , p. 012811 , 2013 . OpenUrl 73. ↵ H. Chen , G. Li , H. Zhang , and Z. Hou , “Optimal allocation of resources for suppressing epidemic spreading on networks ,” Physical Review E , vol. 96 , no. 1 , p. 012321 , 2017 . OpenUrl 74. ↵ M. T. Pearce , A. Agarwala , and D. S. Fisher , “Stabilization of extensive fine-scale diversity by ecologically driven spatiotemporal chaos ,” Proceedings of the National Academy of Sciences , vol. 117 , no. 25 , pp. 14572 – 14583 , 2020 . OpenUrl Abstract / FREE Full Text 75. ↵ J. Hu , Emergent behaviors in complex microbial ecosystems . PhD thesis, Massachusetts Institute of Technology , 2024 . 76. ↵ E. Estrada and J. A. Rodriguez-Velazquez , “Subgraph centrality in complex networks ,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics , vol. 71 , no. 5 , p. 056103 , 2005 . OpenUrl 77. ↵ L. Katz , “A new status index derived from sociometric analysis ,” Psychometrika , vol. 18 , no. 1 , pp. 39 – 43 , 1953 . OpenUrl CrossRef Web of Science 78. ↵ M. E. Newman , “The mathematics of networks ,” The new palgrave encyclopedia of economics , vol. 2 , no. 2008 , pp. 1 – 12 , 2008 . OpenUrl 79. ↵ R. Poulin , M.-C. Boily , and B. R. Mâsse , “Dynamical systems to define centrality in social networks ,” Social networks , vol. 22 , no. 3 , pp. 187 – 220 , 2000 . OpenUrl CrossRef 80. ↵ P. Grindrod and D. J. Higham , “A dynamical systems view of network centrality ,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences , vol. 470 , no. 2165 , p. 20130835 , 2014 . OpenUrl CrossRef 81. ↵ J.-G. Liu , J.-H. Lin , Q. Guo , and T. Zhou , “Locating influential nodes via dynamics-sensitive centrality ,” Scientific reports , vol. 6 , no. 1 , p. 21380 , 2016 . OpenUrl CrossRef PubMed 82. ↵ D.-W. Huang and Z.-G. Yu , “Dynamic-sensitive centrality of nodes in temporal networks ,” Scientific reports , vol. 7 , no. 1 , p. 41454 , 2017 . OpenUrl CrossRef PubMed 83. ↵ A. J. Smola and R. Kondor , “ Kernels and regularization on graphs ,” in Learning Theory and Kernel Machines: 16th Annual Conference on Learning Theory and 7th Kernel Workshop, COLT/Kernel 2003, Washington, DC, USA, August 24-27, 2003. Proceedings , pp. 144 – 158 , Springer , 2003 . 84. ↵ L. Dai , D. Vorselen , K. S. Korolev , and J. Gore , “Generic indicators for loss of resilience before a tipping point leading to population collapse ,” Science , vol. 336 , no. 6085 , pp. 1175 – 1177 , 2012 . OpenUrl Abstract / FREE Full Text 85. J. Karslake , J. Maltas , P. Brumm , and K. B. Wood , “Population density modulates drug inhibition and gives rise to potential bistability of treatment outcomes for bacterial infections ,” PLoS computational biology , vol. 12 , no. 10 , p. e1005098 , 2016 . OpenUrl CrossRef PubMed 86. J. B. Deris , M. Kim , Z. Zhang , H. Okano , R. Hermsen , A. Groisman , and T. Hwa , “The innate growth bistability and fitness landscapes of antibiotic-resistant bacteria ,” Science , vol. 342 , no. 6162 , p. 1237435 , 2013 . OpenUrl Abstract / FREE Full Text 87. N. Frenkel , R. Saar Dover , E. Titon , Y. Shai , and V. Rom-Kedar , “Bistable bacterial growth dynamics in the presence of antimicrobial agents ,” Antibiotics , vol. 10 , no. 1 , p. 87 , 2021 . OpenUrl CrossRef PubMed 88. ↵ E. S. Wright and K. H. Vetsigian , “Inhibitory interactions promote frequent bistability among competing bacteria ,” Nature communications , vol. 7 , no. 1 , p. 11274 , 2016 . OpenUrl CrossRef PubMed 89. ↵ A. Limdi , A. Pérez-Escudero , A. Li , and J. Gore , “Asymmetric migration decreases stability but increases resilience in a heterogeneous metapopulation ,” Nature Communications , vol. 9 , no. 1 , p. 2969 , 2018 . OpenUrl CrossRef PubMed 90. ↵ A. Mishra , P. P. Chakraborty , and S. Dey , “Dispersal evolution diminishes the negative density dependence in dispersal ,” Evolution , vol. 74 , no. 9 , pp. 2149 – 2157 , 2020 . OpenUrl CrossRef PubMed 91. P. P. Chakraborty , L. R. Nemzer , and R. Kassen , “Experimental evidence that network topology can accelerate the spread of beneficial mutations ,” Evolution Letters , vol. 7 , no. 6 , pp. 447 – 456 , 2023 . OpenUrl CrossRef PubMed 92. P. Chakraborty and R. Kassen , “Rapid adaptation and increased genetic parallelism in experimental metapopulations of pseudomonas aeruginosa ,” bioRxiv , pp. 2024 – 10 , 2024 . 93. ↵ P. P. Chakraborty , The Impact of Spatial Network Topologies on the Adaptive Evolution of Pseudomonas aeruginosa Metapopulations . PhD thesis, Université d’Ottawa— University of Ottawa , 2024 . 94. ↵ N. Farrokhian , J. Maltas , M. Dinh , A. Durmaz , P. Ellsworth , M. Hitomi , E. McClure , A. Marusyk , A. Kaznatcheev , and J. G. Scott , “Measuring competitive exclusion in non–small cell lung cancer ,” Science Advances , vol. 8 , no. 26 , p. eabm7212 , 2022 . OpenUrl CrossRef PubMed 95. ↵ J. Maltas , D. S. Tadele , A. Durmaz , C. D. McFarland , M. Hinczewski , and J. G. Scott , “Frequency-dependent ecological interactions increase the prevalence, and shape the distribution, of preexisting drug resistance ,” PRX Life , vol. 2 , no. 2 , p. 023010 , 2024 . OpenUrl CrossRef 96. ↵ P. Newton and Y. Ma , “Maximizing cooperation in the prisoner’s dilemma evolutionary game via optimal control ,” Physical Review E , vol. 103 , no. 1 , p. 012304 , 2021 . OpenUrl 97. R. Dua , Y. Ma , and P. K. Newton , “Are adaptive chemotherapy schedules robust? a three-strategy stochastic evolutionary game theory model ,” Cancers , vol. 13 , no. 12 , p. 2880 , 2021 . OpenUrl CrossRef PubMed 98. K. Stuckey , R. Dua , Y. Ma , J. Parker , and P. Newton , “Optimal dynamic incentive scheduling for hawk-dove evolutionary games ,” Physical Review E , vol. 105 , no. 1 , p. 014412 , 2022 . OpenUrl 99. B. Wölfl , H. Te Rietmole , M. Salvioli , A. Kaznatcheev , F. Thuijsman , J. S. Brown , B. Burgering , and K. StaňkovÁ , “The contribution of evolutionary game theory to understanding and treating cancer ,” Dynamic Games and Applications , vol. 12 , no. 2 , pp. 313 – 342 , 2022 . OpenUrl CrossRef 100. ↵ J. Park and P. Newton , “Stochastic competitive release and adaptive chemotherapy ,” Physical Review E , vol. 108 , no. 3 , p. 034407 , 2023 . OpenUrl View the discussion thread. Back to top Previous Next Posted June 17, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Deciphering and steering population-level response under spatial drug heterogeneity on microhabitat structures Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Deciphering and steering population-level response under spatial drug heterogeneity on microhabitat structures Zhijian Hu , Kevin Wood bioRxiv 2025.02.13.638200; doi: https://doi.org/10.1101/2025.02.13.638200 Share This Article: Copy Citation Tools Deciphering and steering population-level response under spatial drug heterogeneity on microhabitat structures Zhijian Hu , Kevin Wood bioRxiv 2025.02.13.638200; doi: https://doi.org/10.1101/2025.02.13.638200 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Biophysics Subject Areas All Articles Animal Behavior and Cognition (7622) Biochemistry (17648) Bioengineering (13870) Bioinformatics (41880) Biophysics (21423) Cancer Biology (18553) Cell Biology (25458) Clinical Trials (138) Developmental Biology (13364) Ecology (19866) Epidemiology (2067) Evolutionary Biology (24290) Genetics (15589) Genomics (22475) Immunology (17711) Microbiology (40326) Molecular Biology (17145) Neuroscience (88471) Paleontology (666) Pathology (2826) Pharmacology and Toxicology (4815) Physiology (7635) Plant Biology (15114) Scientific Communication and Education (2044) Synthetic Biology (4286) Systems Biology (9815) Zoology (2268)