{"id":"dc759794-4905-4071-b4d3-dd921753ee29","arxiv_id":"2411.16362","paper_version":2,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"A mathematical model predicts that optimal multi-drug therapy switching depends on initial efficacy and pathogen mutation rate, and that constraints on switching create non-trivial optimal protocols.","lead":"This paper builds a stochastic model of how drug efficacy evolves during multi-drug therapy, with therapy switches modeled as random resets, and derives formulas for the average time until resistance develops. It finds that more drugs and faster switching are not always better; with limited or costly switches, there is an optimal protocol that depends on pathogen mutation rate.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Eq. (3)'s parameter d is defined inconsistently between the main text and Appendix B, changing the Bessel order and hence every reported RDT curve and phase diagram.","rationale":"The reader's weakest_assumption focuses on the rotational-symmetry and no-interaction assumptions, which are real limitations but are acknowledged by the authors and are separate from internal consistency. The reader also notes an inconsistency in the printed definition of d, but does not make it the load-bearing issue. My stress-test identifies the d inconsistency as more immediate: even if one accepts the symmetry reduction, Eq. (3) as printed does not specify a unique result because the d used in the Bessel functions is defined two different ways. This is not a disagreement with mainstream biology; it is a checkable internal derivation failure. The paper may well have used a consistent d in its code, in which case the omission is textual, but the printed central claim is not reproducible from the equations as written. Recomputing Fig. 3 with the candidate definitions against the paper's own stochastic simulations settles the question. I therefore keep the reader's CONDITIONAL verdict rather than escalating, because the underlying model and simulations may still support the qualitative conclusions once the definitional issue is resolved. The concrete test also respects the paper's independent support: simulations are used as ground truth, so comparing Eq. (3) against them is the natural adjudicator.","tokens_in":19274,"tokens_out":9937,"duration_ms":102708,"concrete_test":"Recompute Fig. 3a.1 with Eq. (3) using, separately, d = 2[(v/D)+NT] as printed in the main text, d = v/D+NT as printed in Appendix B, and d = NT - v/D as implied by Eq. (B8) with f = -v/eta^2. Overlay these three analytical curves on the Euler-Maruyama simulation of Eq. (2) with v = -8e-5 d^-1, D = 1e-4 d^-1, tau = 10 yr, NT = 2. Whichever d reproduces the simulated mean RDT is the correct implementation; if none matches, Eq. (3) is not the solution of the stated backward equation and the analytical claims in Figs. 3, 5, and 6 need revision.","verdict_should_be":"CONDITIONAL","load_bearing_attack":"The central analytical claim is Eq. (3), the mean RDT for the coupled continuous model. But the paper defines d in two incompatible ways. In the main text, Eq. (3) is introduced with d = 2[(v/D)+NT], and (v/D)+NT is called an effective dimension. In Appendix B, the derivation leading to Eq. (3) states d = v/D+NT and obtains the Bessel equation with index beta = (d/2)-1 (Eqs. B10-B12). Substituting one definition into the Bessel functions J_{d-1}, Y_{d-1} in Eq. (3) gives a different equation than substituting the other, so the printed result is not uniquely specified. The inconsistency is not cosmetic: with the paper's own parameters v = -8e-5 d^-1 and D = 1e-4 d^-1, v/D = -0.8, so the candidate d values differ by roughly a factor of two for NT = 2 and by several units for larger NT. Even the sign convention is suspect: reducing Eq. (B8) with f(eta) = -v/eta^2 gives a radial drift D(NT-1-v/D)/eta, which fixes d = NT - v/D, not d = v/D+NT, unless an unstated sign convention for v is used. Because Eq. (3) is used for Fig. 3, and Eq. (7) built on Eq. (3) determines the Fig. 5 phase diagram, every quantitative optimal-switching claim inherits the ambiguity. The reader's condition is therefore well placed, but the concrete blockage is internal to the derivation, not only an external biological assumption.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","summary":"The manuscript proposes a two-scale stochastic model in which the efficacy of each drug in a multi-drug therapy evolves as a diffusion process in an NT-dimensional efficacy space, therapy switches are modeled as stochastic resets, and drug resistance development is a first-passage event occurring when the radial efficacy sqrt(sum_i eta_i^2) falls below eta_min = 0.4. The authors derive an analytical mean resistance development time (RDT) for a coupled continuous model (Eq. 3) using the backward Fokker-Planck equation and Bessel functions, and a discrete-space expression (Eq. 6) for an uncoupled Markov-chain model. They then study how the mean RDT depends on the therapy switching rate and the number of drugs, identify a threshold initial efficacy separating beneficial from detrimental switching, and examine optimal protocols under therapy constraints such as a minimum switch interval, a limited number of switches, and switching costs.","tokens_in":19704,"tokens_out":6503,"duration_ms":56538,"significance":"If the central formulas are correct, the paper offers a tractable analytical framework for a problem that is usually treated numerically, and it has the notable strengths of checking the continuous formulas against Euler-Maruyama simulations and the discrete model against Gillespie simulations, and of producing falsifiable predictions such as the phase diagram separating 'switch as often as possible' from 'do not switch'. The identification of a threshold initial efficacy separating beneficial from detrimental switching is a useful qualitative insight. However, the main analytical result, Eq. (3), is not uniquely defined because of an inconsistent effective-dimension parameter, so the quantitative predictions, the threshold curves, and the phase diagram need to be revisited before the paper's central claims can be accepted.","major_comments":[{"comment":"The parameter d is defined inconsistently between the main text and the appendix. The main text defines d = 2[(v/D)+NT], while Appendix B3 defines d = v/D+NT and obtains a Bessel equation with index beta = d/2 - 1 (Eq. B12) and a solution involving J_beta, Y_beta, J_{beta+1}, Y_{beta+1} (Eqs. B17-B18). Eq. (3) prints Bessel orders d-1 and d, which do not match the appendix orders under either definition: with the main-text d, the appendix orders would be d/2-1 and d/2, and with the appendix d, the orders would again be d/2-1 and d/2. With v = -8e-5 d^-1 and D = 1e-4 d^-1, the two candidate values of d for NT=2 differ substantially (for example, d = 1.2 or d = 2.4), so Fig. 3, Eq. (7), and the Fig. 5 phase diagram are not uniquely determined. The authors should fix a single definition of d and verify that the Bessel-order reduction in Appendix B reproduces exactly the formula printed as Eq. (3).","section":"Section II.A, Eq. (3), and Appendix B3"},{"comment":"The sign convention for the drift leads to a further discrepancy in d. From Eq. (B8) with f(eta) = -v/eta^2, the radial drift is [D(NT-1) - v]/eta, so the backward operator (B10) has d = NT - v/D if v is taken as a positive magnitude, or d = NT + |v|/D if v is the signed negative value used in the figures. The text instead states d = v/D + NT in Appendix B and d = 2[(v/D)+NT] in the main text. With v = -8e-5 and D = 1e-4, this sign convention changes d by 1.6 for NT=2, which is the same order as the values of d themselves. Please clarify the sign convention for v and re-derive d directly from Eq. (B8).","section":"Appendix B, Eqs. (B7)-(B10), and Eq. (2)"},{"comment":"Several parameters needed to reproduce the simulations are not specified numerically. In particular, the reflecting boundary eta_max used in Eq. (3) and in boundary condition (B15) is never given a value; the number of lattice states M in the uncoupled discrete model is not listed for Figs. 4 and 6; and the cost parameter c in Eq. (9) is not given for the uncoupled-model curves in Fig. 6c. Without these values, the agreement shown between analytics and simulation cannot be checked, and the reader cannot assess whether M is large enough for the Kramers-Moyal expansion of Appendix C1 to be valid.","section":"Figs. 3-6 and Section III"},{"comment":"The absorbing threshold eta_min = 0.4 is introduced as a fixed ad hoc value and is used in both the continuous and discrete models. Because eta_min defines the absorbing boundary, the threshold initial efficacy eta_th, the areas S+ and S- in Fig. 3d, and the phase boundary in Fig. 5c are all functions of this choice. A sensitivity analysis over a plausible range of eta_min, or a calibration of eta_min to the host-pathogen model of Appendix A, is needed to establish that the qualitative phase diagram is robust to this assumption.","section":"Section II, Eq. defining eta_min"}],"minor_comments":[{"comment":"The phrase 'an NT-dimensional lattice of M NT × M NT states' should read 'M^NT states'; the superscripts are missing.","section":"Section IIB"},{"comment":"The word 'hipersphere' should be 'hypersphere'.","section":"Appendix D"},{"comment":"In the sentence introducing Eq. (3), 'Bessel functions of ordern' should read 'Bessel functions of order n'.","section":"Section IIA"},{"comment":"The denominator eta_NT_max - eta_NT_min is missing superscripts and should read eta_max^NT - eta_min^NT.","section":"Eq. (7)"},{"comment":"The statement that successive switches 'diminishes the therapy switching rate 1/tau -> 0' should specify that this is in the limit gamma -> infinity.","section":"Section IVA"},{"comment":"The caption for panel (d) refers to 'a cost function shown in Eq. (10)', but Eq. (10) is introduced for the coupled model with c = 10^{NT-1}; the uncoupled-model curve in panel (c) uses Eq. (9), and the value of c for that curve is not stated.","section":"Fig. 6 caption"}],"recommendation":"major_revision","confidential_remarks":"The manuscript is within scope for a physics-biology interface journal. The central problem is an internal inconsistency in the derivation of Eq. (3), not a disagreement with consensus biological assumptions. The paper should not be rejected on the biological modeling choices alone, but the analytical formula must be made unambiguous and the missing simulation parameters supplied before the results can be verified."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Colleague — the short version: this is a genuinely useful modeling framework for multi-drug switching, but the central analytical result, Eq. (3), has an internal inconsistency in the definition of d that propagates into every quantitative curve and the phase diagram. The paper is worth refereeing, but only after the derivation is cleaned up.\n\nWhat is actually new: the authors extend their earlier single-therapy stochastic resetting model [44] to NT simultaneous drugs, reduce the multi-dimensional efficacy process to a radial coordinate, and derive a Bessel-function expression for mean resistance development time. They also build an uncoupled lattice model solved via the survival matrix, and use a Kramers-Moyal expansion to relate lattice and continuum parameters. On the applied side, they add realistic constraints — minimum switching intervals, finite number of switches, and switching costs — and identify a phase transition between switching as fast as allowed and not switching at all. That phase diagram in Fig. 5 is a real contribution to therapy design.\n\nWhat the paper does well: the two-scale coupling to host-pathogen dynamics is clearly motivated, and the authors are honest about the independence assumption and the lack of empirical validation. The discrete model with Gillespie simulations is a useful complement to the continuum result. Credit is due for framing the problem as maximizing first-passage time rather than minimizing it, which is the right objective here.\n\nWhere the paper falls down: the stress-test note is correct. In the main text Eq. (3) is introduced with d = 2[(v/D)+NT]; in Appendix B3, d = v/D + NT. These differ by roughly a factor of two for the parameters used. Worse, the radial drift in Eq. (B8), with f(η) = -v/η^2, gives D(NT-1-v/D)/η, which fixes d = NT - v/D, not v/D+NT. So the Bessel order printed in Eq. (3) cannot be reproduced from the stated derivation. Since Eq. (3) feeds directly into Figs. 3 and 5, the quantitative claims are not uniquely specified as printed. This is not cosmetic; a referee would need to see the corrected derivation and rerun figures. Also, several simulation parameters are missing: eta_max, M, and the cost parameter c are never given, so the numerics are not reproducible. The threshold eta_min=0.4 is ad hoc, though acknowledged.\n\nWho this is for: mathematical biologists and clinicians interested in antimicrobial resistance and therapy scheduling. It deserves a serious referee, but the review must focus on the derivation of Eq. (3) and on parameter reporting. I would send it to peer review with a request for major revision, not desk reject; the framework is valuable even if the printed equations need fixing.","headline":"Promising multi-drug switching framework, but the central Eq. (3) has an internal inconsistency in the definition of d that must be fixed before the quantitative results can be trusted.","tokens_in":20171,"tokens_out":7987,"would_cite":false,"duration_ms":64217,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["60J70","60J28","92C60"],"pacs":[],"model":"deepseek-v4-flash","headline":"A stochastic two-scale model derives exact mean times to drug resistance and maps when therapy switching helps, when it hurts, and when the best policy is to do nothing.","keywords":["antimicrobial resistance","therapy switching","stochastic resetting","first passage time","multi-drug therapy","resistance development time","master equation","Bessel functions"],"falsifier":"Run the same first-passage calculation with drug-specific diffusion coefficients $D_i$ and drifts $v_i$ while keeping all other assumptions; if the mean RDT deviates from Eq. (3) by more than the Monte Carlo error at any $N_T$ and $\\tau$, the rotational-symmetry reduction is false. Alternatively, fit Eq. (3) to clinical RDT data and check whether the inferred $\\eta_{\\min}$ is stable across regimens, since instability would show the absorbing boundary is not a fixed model parameter.","tokens_in":19119,"feed_emoji":"💊","tokens_out":5382,"duration_ms":48539,"temperature":0.7,"pith_summary":"This paper seeks to establish that the mean time until a multi-drug therapy fails—the resistance development time (RDT)—can be computed analytically for a wide class of stochastic therapy models. The authors model each drug's efficacy as a coordinate in an $N_T$-dimensional space, treat therapy switches as stochastic resets, and define the absorbing boundary as the development of resistance. Their formulas reveal a threshold initial efficacy separating regimens where switching extends the RDT from those where switching shortens it, and a phase transition between switching as often as allowed and never switching. A sympathetic reader would care because these predictions are concrete enough to be tested against clinical and evolutionary data, and because the phase boundary could directly guide protocol design.","feed_headline":"Exact formula finds when drug switching helps or hurts","feed_subtitle":"Stochastic model shows a threshold efficacy and a phase transition between switching constantly and never switching.","key_machinery":"The load-bearing object is the reduction of an $N_T$-dimensional isotropic diffusion with stochastic resetting to the radial coordinate $\\eta$, exploiting rotational symmetry and a radial drift $f(\\eta) = -v/\\eta^2$. This reduction converts the backward Fokker-Planck equation into a Bessel equation whose solution, Eq. (3), gives the conditioned mean RDT. The discrete model instead builds a transition matrix $W$ on an $N_T$-dimensional lattice and reads the mean absorption time from the survival matrix $S$ in Eq. (6). The Bessel solution and the survival-matrix formula carry every later conclusion about thresholds, phases, and optimal protocols.","core_discovery":"The paper's central claim is that the RDT in multi-drug therapy is captured by two analytically tractable models: a coupled continuous model in which the radial efficacy $\\eta = \\sqrt{\\sum_i \\eta_i^2}$ follows a Bessel-function mean first-passage formula (Eq. 3), and an uncoupled discrete model in which the master equation on an $M^{N_T}$ lattice gives the mean RDT through a survival-matrix formula (Eq. 6). Using these expressions, the authors identify a threshold initial efficacy $\\eta_{\\mathrm{th}}$ that separates parameter regions where therapy switches lengthen the RDT from regions where they shorten it. They also identify a phase boundary in the $(N_T,\\tau_{\\min})$ plane between the strategies \"switch as fast as allowed\" and \"never switch.\" Under limited or costly switching, the optimal switching rate is non-monotonic and depends on the pathogen's mutation rate.","pith_inferences":["Editorial extension: the rotational-symmetry reduction suggests a testable collapse—for a fixed drug count, measured RDTs at different initial efficacies should fall on the same Bessel curve after rescaling by $\\lambda = \\sqrt{D\\tau}$, and deviations would expose drug-specific or interaction effects.","Editorial extension: the phase diagram implies a cheap clinical prescription: when a regimen already uses many drugs, adding more switching may be wasted effort, so determining which side of the $(N_T,\\tau_{\\min})$ boundary a regimen lies on could guide protocol design without new simulations.","Editorial extension: the model's $\\eta_{\\min}=0.4$ is an ad hoc threshold; calibrating the failure boundary from within-host or clinical data would convert the predicted phase boundary into a quantitative, testable prediction with error bars."],"forward_implications":["For patients whose initial efficacy is above $\\eta_{\\mathrm{th}}$, increasing the switching rate or the number of drugs extends the mean RDT; below that threshold, switching shortens the RDT, so therapy switching can be actively harmful.","With a fixed minimum interval between switches, the optimal policy is either to switch as often as allowed or never, and the choice is decided by a phase boundary in the drug-count versus minimum-interval plane.","When only a limited number of switches is available, the mean RDT is non-monotonic in the switching rate, so an optimal finite switching rate exists.","Pathogens with a larger mutation rate (diffusion constant $D$) benefit from more simultaneous drugs, while slowly mutating pathogens do better with fewer drugs and fewer switches.","Because the never-switch phase also avoids therapy costs, increasing the number of simultaneous drugs can simultaneously maximize the RDT and reduce treatment expense."],"supporting_citations":[{"why":"Introduces stochastic resetting, the mechanism by which therapy switches return efficacy to its initial value; the paper's dynamics are built on this framework.","marker":"[31]"},{"why":"Previous single-drug stochastic-resetting model of therapy administration that this work extends to $N_T$ simultaneous drugs; supplies the baseline result and parameters.","marker":"[44]"},{"why":"Within-host HIV-1 infection model whose steady state fixes the infection-rate coupling and the absorbing/reflecting boundary geometry of efficacy space.","marker":"[23]"},{"why":"Provides the backward Fokker-Planck formalism used to derive the Bessel-function mean first-passage formula in Eq. (3).","marker":"[46]"},{"why":"Gillespie algorithm used to simulate the uncoupled discrete model and to verify the analytical mean RDT and optimal-switching results.","marker":"[49]"},{"why":"Supplies the survival-matrix method used in Eq. (6) to compute mean residence times before absorption on the lattice.","marker":"[51–53]"}],"fun_headline_variants":["Drug-switch threshold: when to switch, when to hold","Phase transition found in drug-switching strategy","Optimal drug-switching rate hinges on a critical efficacy","Exact formulas reveal when switching drugs helps most"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The exact results assume every drug acts identically and independently: isotropic, drug-agnostic diffusion, a rotationally symmetric drift toward failure, and an absorbing boundary fixed at the hand-chosen value $\\eta_{\\min}=0.4$; break that symmetry or change the threshold, and the formulas, thresholds, and phase diagram all shift.","fun_headline_variants_meta":{"raw":{"variants":["Drug-switch threshold: when to switch, when to hold","Phase transition found in drug-switching strategy","Optimal drug-switching rate hinges on a critical efficacy","Exact formulas reveal when switching drugs helps most"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000277,"raw_usage":{"total_tokens":1655,"prompt_tokens":953,"completion_tokens":702,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":569,"completion_tokens_details":{"reasoning_tokens":639}},"tokens_in":569,"tokens_out":702,"duration_ms":7642,"temperature":1.0,"reasoning_tokens":639,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-12T13:12:17.344092+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"Run the same first-passage calculation with drug-specific diffusion coefficients $D_i$ and drifts $v_i$ while keeping all other assumptions; if the mean RDT deviates from Eq. (3) by more than the Monte Carlo error at any $N_T$ and $\\tau$, the rotational-symmetry reduction is false. Alternatively, fit Eq. (3) to clinical RDT data and check whether the inferred $\\eta_{\\min}$ is stable across regimens, since instability would show the absorbing boundary is not a fixed model parameter.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces stochastic resetting, the mechanism by which therapy switches return efficacy to its initial value; the paper's dynamics are built on this framework."},{"cited_title":"Ramoso, J","cited_arxiv_id":null,"evidence_quote":"Previous single-drug stochastic-resetting model of therapy administration that this work extends to $N_T$ simultaneous drugs; supplies the baseline result and parameters."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Within-host HIV-1 infection model whose steady state fixes the infection-rate coupling and the absorbing/reflecting boundary geometry of efficacy space."},{"cited_title":"Gardiner, in Stochastic Methods, Springer series in synergetics (Springer Berlin Heidelberg, Berlin, Heidel- berg, 2008) pp","cited_arxiv_id":null,"evidence_quote":"Provides the backward Fokker-Planck formalism used to derive the Bessel-function mean first-passage formula in Eq. (3)."}],"review_version":1}