Novel optimization techniques for inferring heterogeneous population dynamics

Chenyu Wu; Nuozhou Wang; Casey Garner; Kevin Leder; Shuzhong Zhang; Chenyu Wu; Nuozhou Wang; Casey Garner; Kevin Leder; Shuzhong Zhang

doi:10.3934/mbe.2026074

Mathematical Biosciences and Engineering

2026, Volume 23, Issue 7: 2018-2054. doi: 10.3934/mbe.2026074

Previous Article Next Article

Research article

Novel optimization techniques for inferring heterogeneous population dynamics

1.
Department of Industrial & Systems Engineering, University of Minnesota, USA
2.
School of Mathematics, University of Minnesota, USA

Received: 20 November 2025 Revised: 19 January 2026 Accepted: 18 March 2026 Published: 09 July 2026

In this paper, we introduce a new optimization algorithm that is well suited to solve parameter estimation problems that arise when inferring heterogeneous population dynamics. In these estimation problems, parameter estimation is complicated by the presence of two types of constraints: inequality constraints (e.g., non-negativity and boundedness of rates (so-called box-constraints)) and equality constraints that arise due to the need of the population fractions to sum to one. We call our new method cubic regularized Newton with affine scaling (CRNAS). In contrast to so-called first-order methods, which solely rely on the gradient of the objective function, our method utilizes the Hessian of the objective. As a result, it is able to focus on points that satisfy the second-order optimality conditions, as opposed to first-order methods that simply converge to critical points. This is an important feature in parameter estimation problems, where the objective function is often non-convex; as a result, there can be many critical points, which makes it nearly impossible to identify the global minimum. We use an affine scaling approach to handle a wide class of constraints, including equality constraints. We establish that CRNAS identifies a point that satisfies $ \epsilon $-approximate second-order optimality conditions within $ O(\epsilon^{-3/2}) $ iterations. Finally, we compare CRNAS with MATLAB's optimization solver fmincon on three different test problems. These test problems all feature mixtures of heterogeneous populations, a problem setting that CRNAS is particularly well-suited for. Our numerical simulations show that CRNAS has a favorable performance, thereby performing comparable, if not better than, fmincon in accuracy and computational cost for most of our examples.
- nonlinear programming,
- nonconvex optimization,
- parameter estimation,
- heterogeneous population dynamics,
- branching process
Citation: Chenyu Wu, Nuozhou Wang, Casey Garner, Kevin Leder, Shuzhong Zhang. Novel optimization techniques for inferring heterogeneous population dynamics[J]. Mathematical Biosciences and Engineering, 2026, 23(7): 2018-2054. doi: 10.3934/mbe.2026074

Related Papers:

Abstract

In this paper, we introduce a new optimization algorithm that is well suited to solve parameter estimation problems that arise when inferring heterogeneous population dynamics. In these estimation problems, parameter estimation is complicated by the presence of two types of constraints: inequality constraints (e.g., non-negativity and boundedness of rates (so-called box-constraints)) and equality constraints that arise due to the need of the population fractions to sum to one. We call our new method cubic regularized Newton with affine scaling (CRNAS). In contrast to so-called first-order methods, which solely rely on the gradient of the objective function, our method utilizes the Hessian of the objective. As a result, it is able to focus on points that satisfy the second-order optimality conditions, as opposed to first-order methods that simply converge to critical points. This is an important feature in parameter estimation problems, where the objective function is often non-convex; as a result, there can be many critical points, which makes it nearly impossible to identify the global minimum. We use an affine scaling approach to handle a wide class of constraints, including equality constraints. We establish that CRNAS identifies a point that satisfies $ \epsilon $-approximate second-order optimality conditions within $ O(\epsilon^{-3/2}) $ iterations. Finally, we compare CRNAS with MATLAB's optimization solver fmincon on three different test problems. These test problems all feature mixtures of heterogeneous populations, a problem setting that CRNAS is particularly well-suited for. Our numerical simulations show that CRNAS has a favorable performance, thereby performing comparable, if not better than, fmincon in accuracy and computational cost for most of our examples.

References

[1]	H. Hass, C. Loos, E. Raimúndez-Álvarez, J. Timmer, J. Hasenauer, C. Kreutz, Benchmark problems for dynamic modeling of intracellular processes, Bioinformatics, 35 (2019), 3073–3082. https://doi.org/10.1093/bioinformatics/btz020 doi: 10.1093/bioinformatics/btz020
[2]	Y. Cao, S. Li, L. Petzold, R. Serban, Adjoint sensitivity analysis for differential-algebraic equations: The adjoint DAE system and its numerical solution, SIAM J. Sci. Comput., 24 (2003), 1076–1089. https://doi.org/10.1137/S1064827501380630 doi: 10.1137/S1064827501380630
[3]	E. G. Birgin, J. M. Martínez, On regularization and active-set methods with complexity for constrained optimization, SIAM J. Optim., 28 (2018), 1367–1395. https://doi.org/10.1137/17M1127107 doi: 10.1137/17M1127107
[4]	E. G. Birgin, G. Haeser, A. Ramos, Augmented Lagrangians with constrained subproblems and convergence to second-order stationary points, Comput. Optim. Appl., 69 (2018), 51–75. https://doi.org/10.1007/s10589-017-9937-2 doi: 10.1007/s10589-017-9937-2
[5]	A. Daigeler, L. Klein-Hitpass, A. M. Chromik, O. Müller, J. Hauser, H. H. Homann, et al., Heterogeneous in vitro effects of doxorubicin on gene expression in primary human liposarcoma cultures, BMC Cancer, 8 (2008), 1–17. https://doi.org/10.1186/1471-2407-8-313 doi: 10.1186/1471-2407-8-313
[6]	M. Nouiehed, J. D. Lee, M. Razaviyayn, Convergence to second-order stationarity for constrained non-convex optimization, preprint, arXiv: 1810.02024. https://doi.org/10.48550/arXiv.1810.02024
[7]	C. Wu, E. B. Gunnarsson, E. M. Myklebust, A. Köhn-Luque, D. S. Tadele, J. M. Enserink, et al., Using birth-death processes to infer tumor subpopulation structure from live-cell imaging drug screening data, PLoS Comput. Biol., 20 (2024), e1011888. https://doi.org/10.1371/journal.pcbi.1011888 doi: 10.1371/journal.pcbi.1011888
[8]	L. Ljung, T. Chen, Convexity issues in system identification, 2013 10th IEEE Int. Conf. Control Autom. (ICCA), (2013), 1–9. https://doi.org/10.1109/ICCA.2013.6565206
[9]	C. G. Moles, P. Mendes, J. R. Banga, Parameter estimation in biochemical pathways: A comparison of global optimization methods, Genome Res., 13 (2003), 2467–2474. https://doi.org/10.1101/gr.1262503 doi: 10.1101/gr.1262503
[10]	L. Schmiester, D. Weindl, J. Hasenauer, Efficient gradient-based parameter estimation for dynamic models using qualitative data, Bioinformatics, 37 (2021), 4493–4500. https://doi.org/10.1093/bioinformatics/btab512 doi: 10.1093/bioinformatics/btab512
[11]	A. Gábor, J. R. Banga, Robust and efficient parameter estimation in dynamic models of biological systems, BMC Syst. Biol., 9 (2015), 1–25. https://doi.org/10.1186/s12918-015-0219-2 doi: 10.1186/s12918-015-0219-2
[12]	F. Fröhlich, P. K. Sorger, Fides: Reliable trust-region optimization for parameter estimation of ordinary differential equation models, PLoS Comput. Biol., 18 (2022), e1010322. https://doi.org/10.1371/journal.pcbi.1010322 doi: 10.1371/journal.pcbi.1010322
[13]	P. Stapor, F. Fröhlich, J. Hasenauer, Optimization and profile calculation of ode models using second order adjoint sensitivity analysis, Bioinformatics, 34 (2018), i151–i159. https://doi.org/10.1093/bioinformatics/bty230 doi: 10.1093/bioinformatics/bty230
[14]	N. Boumal, V. Voroninski, A. Bandeira, The non-convex Burer-Monteiro approach works on smooth semidefinite programs, in 30th Conference on Neural Information Processing Systems (NIPS 2016), (2016).
[15]	R. Ge, J. D. Lee, T. Ma, Matrix completion has no spurious local minimum, in Advances in Neural Information Processing Systems 29 (NIPS 2016), (2016).
[16]	F. E. Curtis, D. P. Robinson, M. Samadi, An inexact regularized Newton framework with a worst-case iteration complexity of for nonconvex optimization, IMA J. Numer. Anal., 39 (2019), 1296–1327. https://doi.org/10.1093/imanum/dry022 doi: 10.1093/imanum/dry022
[17]	M. O'Neill, S. J. Wright, A log-barrier Newton-CG method for bound constrained optimization with complexity guarantees, IMA J. Numer. Anal., 41 (2021), 84–121. https://doi.org/10.1093/imanum/drz074 doi: 10.1093/imanum/drz074
[18]	C. W. Royer, M. O'Neill, S. J. Wright, A Newton-CG algorithm with complexity guarantees for smooth unconstrained optimization, Math. Program., 180 (2020), 451–488. https://doi.org/10.1007/s10107-019-01362-7 doi: 10.1007/s10107-019-01362-7
[19]	P. Dvurechensky, M. Staudigl, Barrier algorithms for constrained non-convex optimization, preprint, arXiv: 2404.18724. https://doi.org/10.48550/arXiv.2404.18724
[20]	C. He, Z. Lu, A Newton-CG based barrier method for finding a second-order stationary point of nonconvex conic optimization with complexity guarantees, SIAM J. Optim., 33 (2023), 1191–1222. https://doi.org/10.1137/21m1457011 doi: 10.1137/21m1457011
[21]	Y. Nesterov, A. Nemirovskii, Interior-point Polynomial Algorithms in Convex Programming, SIAM, 1994. https://doi.org/10.1137/1.9781611970791
[22]	Y. Nesterov, B. T. Polyak, Cubic regularization of Newton method and its global performance, Math. Program., 108 (2006), 177–205. https://doi.org/10.1007/s10107-006-0706-8 doi: 10.1007/s10107-006-0706-8
[23]	A. Griewank, The modification of Newton's method for unconstrained optimization by bounding cubic terms, Technical Report NA/12, (1981). https://doi.org/10.13140/RG.2.1.4097.2960
[24]	J. Lagarias, R. Vanderbei, Ⅱ Dikin's convergence result for the affine scaling algorithm, Contemp. Math, 114 (1990), 109–119. https://doi.org/10.1090/conm/114/1097868 doi: 10.1090/conm/114/1097868
[25]	E. R. Barnes, A variation on Karmarkar's algorithm for solving linear programming problems, Math. Program., 36 (1986), 174–182. https://doi.org/10.1007/bf02592024 doi: 10.1007/bf02592024
[26]	R. J. Vanderbei, M. S. Meketon, B. A. Freedman, A modification of Karmarkar's linear programming algorithm, Algorithmica, 1 (1986), 395–407. https://doi.org/10.1007/BF01840454 doi: 10.1007/BF01840454
[27]	MATLAB Optimization Toolbox, MATLAB Optimization Toolbox, R2022a. The MathWorks, Natick, MA, USA. Available from: https://www.mathworks.com.
[28]	E. B. Gunnarsson, J. Foo, K. Leder, Statistical inference of the rates of cell proliferation and phenotypic switching in cancer, J. Theor. Biol., 568 (2023), 111497. https://doi.org/10.1016/j.jtbi.2023.111497 doi: 10.1016/j.jtbi.2023.111497
[29]	A. F. Villaverde, F. Fröhlich, D. Weindl, J. Hasenauer, J. R. Banga, Benchmarking optimization methods for parameter estimation in large kinetic models, Bioinformatics, 35 (2019), 830–838. https://doi.org/10.1093/bioinformatics/bty736 doi: 10.1093/bioinformatics/bty736
[30]	D. K. Agrawal, B. J. Smith, P. D. Sottile, G. Hripcsak, D. J. Albers, Quantifiable identification of flow-limited ventilator dyssynchrony with the deformed lung ventilator model, Comput. Biol. Med., 173 (2024), 108349. https://doi.org/10.1016/j.compbiomed.2024.108349 doi: 10.1016/j.compbiomed.2024.108349
[31]	A. Köhn-Luque, E. M. Myklebust, D. S. Tadele, M. Giliberto, L. Schmiester, J. Noory, et al., Phenotypic deconvolution in heterogeneous cancer cell populations using drug-screening data, Cell Rep. Methods, 3 (2023), 100417. https://doi.org/10.1016/j.crmeth.2023.100417 doi: 10.1016/j.crmeth.2023.100417
[32]	R. Gesztelyi, J. Zsuga, A. Kemeny-Beke, B. Varga, B. Juhasz, A. Tosaki, The Hill equation and the origin of quantitative pharmacology, Arch. Hist. Exact Sci., 66 (2012), 427–438. https://doi.org/10.1007/s00407-012-0098-5 doi: 10.1007/s00407-012-0098-5
[33]	A. V. Hill, The possible effects of the aggregation of the molecules of hemoglobin on its dissociation curves, J. Physiol., 40 (1910), iv–vii.
[34]	S. R. Gadagkar, G. B. Call, Computational tools for fitting the Hill equation to dose–response curves, J. Pharmacol. Toxicol. Methods, 71 (2015), 68–76. https://doi.org/10.1016/j.vascn.2014.08.006 doi: 10.1016/j.vascn.2014.08.006
[35]	W. Jin, S. W. McCue, M. J. Simpson, Extended logistic growth model for heterogeneous populations, J. Theor. Biol., 445 (2018), 51–61. https://doi.org/10.1016/j.jtbi.2018.02.027 doi: 10.1016/j.jtbi.2018.02.027
[36]	C. Cartis, N. I. M. Gould, P. L. Toint, Adaptive cubic regularisation methods for unconstrained optimization. Part Ⅰ: Motivation, convergence and numerical results, Math. Program., 127 (2011), 245–295. https://doi.org/10.1007/s10107-009-0286-5 doi: 10.1007/s10107-009-0286-5
[37]	A. R. Conn, N. I. M. Gould, P. L. Toint, Trust region methods, Soc. Ind. Appl. Math., (2000). https://doi.org/10.1137/1.9780898719857
[38]	D. T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, J. Comput. Phys., 22 (1976), 403–434. https://doi.org/10.1016/0021-9991(76)90041-3 doi: 10.1016/0021-9991(76)90041-3

Reader Comments

Your name:*

Email:*
© 2026 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)