Coarse-graining Hamiltonian systems using WSINDy.

Weak form equation learning and surrogate modeling has proven to be computationally efficient and robust to measurement noise in a wide range of applications including ODE, PDE, and SDE discovery, as well as in coarse-graining applications, such as homogenization and mean-field descriptions of interacting particle systems. In this work we extend this coarse-graining capability to the setting of Hamiltonian dynamics which possess approximate symmetries associated with timescale separation. A smooth ε -dependent Hamiltonian vector field X ε possesses an approximate symmetry if the limiting vector field X 0 = lim ε → 0 X ε possesses an exact symmetry. Such approximate symmetries often lead to the existence of a Hamiltonian system of reduced dimension that may be used to efficiently capture the dynamics of the symmetry-invariant dependent variables. Deriving such reduced systems, or approximating them numerically, is an ongoing challenge. We demonstrate that WSINDy can successfully identify this reduced Hamiltonian system in the presence of large perturbations imparted in the ε > 0 regime, while remaining robust to extrinsic noise. This is significant in part due to the nontrivial means by which such systems are derived analytically. WSINDy naturally preserves the Hamiltonian structure by restricting to a trial basis of Hamiltonian vector fields. The methodology is computationally efficient, often requiring only a single trajectory to learn the global reduced Hamiltonian, and avoiding forward solves in the learning process. In this way, we argue that weak-form equation learning is particularly well-suited for Hamiltonian coarse-graining. Using nearly-periodic Hamiltonian systems as a prototypical class of systems with approximate symmetries, we show that WSINDy robustly identifies the correct leading-order system, with dimension reduced by at least two, upon observation of the relevant degrees of freedom. While our main contribution is computational, we also provide a contribution to the literature on averaging theory by proving that first-order averaging at the level of vector fields preserves Hamiltonian structure in nearly-periodic Hamiltonian systems. This provides theoretical justification for our approach as WSINDy's computations occur at the level of Hamiltonian vector fields. We illustrate the efficacy of our proposed method using physically relevant examples, including coupled oscillator dynamics, the Hénon-Heiles system for stellar motion within a galaxy, and the dynamics of charged particles.

Direct Estimation of Parameters in ODE Models Using WENDy: Weak-Form Estimation of Nonlinear Dynamics.

We introduce the Weak-form Estimation of Nonlinear Dynamics (WENDy) method for estimating model parameters for non-linear systems of ODEs. Without relying on any numerical differential equation solvers, WENDy computes accurate estimates and is robust to large (biologically relevant) levels of measurement noise. For low dimensional systems with modest amounts of data, WENDy is competitive with conventional forward solver-based nonlinear least squares methods in terms of speed and accuracy. For both higher dimensional systems and stiff systems, WENDy is typically both faster (often by orders of magnitude) and more accurate than forward solver-based approaches. The core mathematical idea involves an efficient conversion of the strong form representation of a model to its weak form, and then solving a regression problem to perform parameter inference. The core statistical idea rests on the Errors-In-Variables framework, which necessitates the use of the iteratively reweighted least squares algorithm. Further improvements are obtained by using orthonormal test functions, created from a set of [Formula: see text] bump functions of varying support sizes.We demonstrate the high robustness and computational efficiency by applying WENDy to estimate parameters in some common models from population biology, neuroscience, and biochemistry, including logistic growth, Lotka-Volterra, FitzHugh-Nagumo, Hindmarsh-Rose, and a Protein Transduction Benchmark model. Software and code for reproducing the examples is available at https://github.com/MathBioCU/WENDy .

Medial Prefrontal Cortex to Medial Septum Pathway Activation Improves Cognitive Flexibility in Rats.

BACKGROUND: The medial prefrontal cortex (mPFC) is necessary for cognitive flexibility and projects to medial septum (MS). MS activation improves strategy switching, a common measure of cognitive flexibility, likely via its ability to regulate midbrain dopamine (DA) neuron population activity. We hypothesized that the mPFC to MS pathway (mPFC-MS) may be the mechanism by which the MS regulates strategy switching and DA neuron population activity. METHODS: Male and female rats learned a complex discrimination strategy across 2 different training time points: a constant length (10 days) and a variable length that coincided with each rat meeting an acquisition-level performance threshold (males: 5.3 ± 0.3 days, females: 3.8 ± 0.3 days). We then chemogenetically activated or inhibited the mPFC-MS pathway and measured each rat's ability to inhibit the prior learned discrimination strategy and switch to a prior ignored discrimination strategy (strategy switching). RESULTS: Activation of the mPFC-MS pathway improved strategy switching after 10 days of training in both sexes. Inhibition of the pathway produced a modest improvement in strategy switching that was quantitatively and qualitatively different from pathway activation. Neither activation nor inhibition of the mPFC-MS pathway affected strategy switching following the acquisition-level performance threshold training regimen. Activation, but not inhibition, of the mPFC-MS pathway bidirectionally regulated DA neuron activity in the ventral tegmental area and substantia nigra pars compacta, similar to general MS activation. CONCLUSIONS: This study presents a potential top-down circuit from the prefrontal cortex to the midbrain by which DA activity can be manipulated to promote cognitive flexibility.

Direct Estimation of Parameters in ODE Models Using WENDy: Weak-form Estimation of Nonlinear Dynamics.

We introduce the Weak-form Estimation of Nonlinear Dynamics (WENDy) method for estimating model parameters for non-linear systems of ODEs. Without relying on any numerical differential equation solvers, WENDy computes accurate estimates and is robust to large (biologically relevant) levels of measurement noise. For low dimensional systems with modest amounts of data, WENDy is competitive with conventional forward solver-based nonlinear least squares methods in terms of speed and accuracy. For both higher dimensional systems and stiff systems, WENDy is typically both faster (often by orders of magnitude) and more accurate than forward solver-based approaches. The core mathematical idea involves an efficient conversion of the strong form representation of a model to its weak form, and then solving a regression problem to perform parameter inference. The core statistical idea rests on the Errors-In-Variables framework, which necessitates the use of the iteratively reweighted least squares algorithm. Further improvements are obtained by using orthonormal test functions, created from a set of C-infinity bump functions of varying support sizes. We demonstrate the high robustness and computational efficiency by applying WENDy to estimate parameters in some common models from population biology, neuroscience, and biochemistry, including logistic growth, Lotka-Volterra, FitzHugh-Nagumo, Hindmarsh-Rose, and a Protein Transduction Benchmark model. Software and code for reproducing the examples is available at (https://github.com/MathBioCU/WENDy).

Medial septum activation improves strategy switching once strategies are well-learned via bidirectional regulation of dopamine neuron population activity.

Strategy switching is a form of cognitive flexibility that requires inhibiting a previously successful strategy and switching to a new strategy of a different categorical modality. It is dependent on dopamine (DA) receptor activation and release in ventral striatum and prefrontal cortex, two primary targets of ventral tegmental area (VTA) DA projections. Although the circuitry that underlies strategy switching early in learning has been studied, few studies have examined it after extended discrimination training. This may be important as DA activity and release patterns change across learning, with several studies demonstrating a critical role for substantia nigra pars compacta (SNc) DA activity and release once behaviors are well-learned. We have demonstrated that medial septum (MS) activation simultaneously increased VTA and decreased SNc DA population activity, as well as improved reversal learning via these actions on DA activity. We hypothesized that MS activation would improve strategy switching both early in learning and after extended training through its ability to increase VTA DA population activity and decrease SNc DA population activity, respectively. We chemogenetically activated the MS of male and female rats and measured their performance on an operant-based strategy switching task following 1, 10, or 15 days of discrimination training. Contrary to our hypothesis, MS activation did not affect strategy switching after 1 day of discrimination training. MS activation improved strategy switching after 10 days of training, but only in females. MS activation improved strategy switching in both sexes after 15 days of training. Infusion of bicuculline into the ventral subiculum (vSub) inhibited the MS-mediated decrease in SNc DA population activity and attenuated the improvement in strategy switching. Intra-vSub infusion of scopolamine inhibited the MS-mediated increase in VTA DA population activity but did not affect the improvement in strategy switching. Intra-vSub infusion of both bicuculline and scopolamine inhibited the MS-mediated effects on DA population activity in both the SNc and VTA and completely prevented the improvement in strategy switching. These data indicate that MS activation improves strategy switching once the original strategy has been sufficiently well-learned, and that this may occur via the MS's regulation of DA neuron responsivity.

Learning mean-field equations from particle data using WSINDy.

We develop a weak-form sparse identification method for interacting particle systems (IPS) with the primary goals of reducing computational complexity for large particle number N and offering robustness to either intrinsic or extrinsic noise. In particular, we use concepts from mean-field theory of IPS in combination with the weak-form sparse identification of nonlinear dynamics algorithm (WSINDy) to provide a fast and reliable system identification scheme for recovering the governing stochastic differential equations for an IPS when the number of particles per experiment N is on the order of several thousands and the number of experiments M is less than 100. This is in contrast to existing work showing that system identification for N less than 100 and M on the order of several thousand is feasible using strong-form methods. We prove that under some standard regularity assumptions the scheme converges with rate O(N-1∕2) in the ordinary least squares setting and we demonstrate the convergence rate numerically on several systems in one and two spatial dimensions. Our examples include a canonical problem from homogenization theory (as a first step towards learning coarse-grained models), the dynamics of an attractive-repulsive swarm, and the IPS description of the parabolic-elliptic Keller-Segel model for chemotaxis. Code is available at https://github.com/MathBioCU/WSINDy_IPS.

ANALYTICAL SINGULAR VALUE DECOMPOSITION FOR A CLASS OF STOICHIOMETRY MATRICES.

We present the analytical singular value decomposition of the stoichiometry matrix for a spatially discrete reaction-diffusion system. The motivation for this work is to develop a matrix decomposition that can reveal hidden spatial flux patterns of chemical reactions. We consider a 1D domain with two subregions sharing a single common boundary. Each of the subregions is further partitioned into a finite number of compartments. Chemical reactions can occur within a compartment, whereas diffusion is represented as movement between adjacent compartments. Inspired by biology, we study both (1) the case where the reactions on each side of the boundary are different and only certain species diffuse across the boundary and (2) the case where reactions and diffusion are spatially homogeneous. We write the stoichiometry matrix for these two classes of systems using a Kronecker product formulation. For the first scenario, we apply linear perturbation theory to derive an approximate singular value decomposition in the limit as diffusion becomes much faster than reactions. For the second scenario, we derive an exact analytical singular value decomposition for all relative diffusion and reaction time scales. By writing the stoichiometry matrix using Kronecker products, we show that the singular vectors and values can also be written concisely using Kronecker products. Ultimately, we find that the singular value decomposition of the reaction-diffusion stoichiometry matrix depends on the singular value decompositions of smaller matrices. These smaller matrices represent modified versions of the reaction-only stoichiometry matrices and the analytically known diffusion-only stoichiometry matrix. Lastly, we present the singular value decomposition of the model for the Calvin cycle in cyanobacteria and demonstrate the accuracy of our formulation. The MATLAB code, available at www.github.com/MathBioCU/ReacDiffStoicSVD, provides routines for efficiently calculating the SVD for a given reaction network on a 1D spatial domain.

Online Weak-form Sparse Identification of Partial Differential Equations.

This paper presents an online algorithm for identification of partial differential equations (PDEs) based on the weak-form sparse identification of nonlinear dynamics algorithm (WSINDy). The algorithm is online in the sense that if performs the identification task by processing solution snapshots that arrive sequentially. The core of the method combines a weak-form discretization of candidate PDEs with an online proximal gradient descent approach to the sparse regression problem. In particular, we do not regularize the ℓ0-pseudo-norm, instead finding that directly applying its proximal operator (which corresponds to a hard thresholding) leads to efficient online system identification from noisy data. We demonstrate the success of the method on the Kuramoto-Sivashinsky equation, the nonlinear wave equation with time-varying wavespeed, and the linear wave equation, in one, two, and three spatial dimensions, respectively. In particular, our examples show that the method is capable of identifying and tracking systems with coefficients that vary abruptly in time, and offers a streaming alternative to problems in higher dimensions.

WEAK SINDY FOR PARTIAL DIFFERENTIAL EQUATIONS.

Sparse Identification of Nonlinear Dynamics (SINDy) is a method of system discovery that has been shown to successfully recover governing dynamical systems from data [6, 39]. Recently, several groups have independently discovered that the weak formulation provides orders of magnitude better robustness to noise. Here we extend our Weak SINDy (WSINDy) framework introduced in [28] to the setting of partial differential equations (PDEs). The elimination of pointwise derivative approximations via the weak form enables effective machine-precision recovery of model coefficients from noise-free data (i.e. below the tolerance of the simulation scheme) as well as robust identification of PDEs in the large noise regime (with signal-to-noise ratio approaching one in many well-known cases). This is accomplished by discretizing a convolutional weak form of the PDE and exploiting separability of test functions for efficient model identification using the Fast Fourier Transform. The resulting WSINDy algorithm for PDEs has a worst-case computational complexity of O ( N D + 1 log ( N ) ) for datasets with N points in each of D + 1 dimensions. Furthermore, our Fourier-based implementation reveals a connection between robustness to noise and the spectra of test functions, which we utilize in an a priori selection algorithm for test functions. Finally, we introduce a learning algorithm for the threshold in sequential-thresholding least-squares (STLS) that enables model identification from large libraries, and we utilize scale invariance at the continuum level to identify PDEs from poorly-scaled datasets. We demonstrate WSINDy's robustness, speed and accuracy on several challenging PDEs. Code is publicly available on GitHub at https://github.com/MathBioCU/WSINDy_PDE.

Protocol for Analysis and Consolidation of TrackMate Outputs for Measuring Two-Dimensional Cell Motility using Nuclear Tracking.

Collective cellular migration plays a key role in many fundamental biological processes including development, wound healing, and cancer metastasis. To understand the regulation of cell motility, we must be able to measure it easily and consistently under different conditions. Here we describe a method for measuring and quantifying single-cell and bulk motility of HaCaT keratinocytes using a nuclear stain. This method includes a MATLAB script for analyzing TrackMate output files to calculate displacements, motility rates, and trajectory angles in single cells and in bulk for an imaging site. This motility analysis script allows for quick, straightforward, and scalable analysis of cell motility rates from TrackMate data and could be broadly used to identify and study the regulation of motility in epithelial cells. We also provide a MATLAB script for reorganizing microscopy videos collected on a microscope and converting them to TIF stacks, which can be analyzed using the ImageJ TrackMate plugin in bulk. Using this methodology to explore the roles of adherens junctions and actin cytoskeletal dynamics in regulating cell motility in HaCaT keratinocytes, we demonstrate evidence that Arp2/3 activity is required for the elevated motility seen after α-catenin depletion in HaCaT keratinocytes.

BOUNDEDNESS OF A CLASS OF SPATIALLY DISCRETE REACTION-DIFFUSION SYSTEMS.

Although the spatially discrete reaction-diffusion equation is often used to describe biological processes, the effect of diffusion in this framework is not fully understood. In the spatially continuous case, the incorporation of diffusion can cause blow-up with respect to the L∞ norm, and criteria exist to determine whether the system is bounded for all time. However, no equivalent criteria exist for the discrete reaction-diffusion system. Due to the possible dynamical differences between these two system types and the advantage of using the spatially discrete representation to describe biological processes, it is worth examining the discrete system independently of the continuous system. Therefore, the focus of this paper is on determining sufficient conditions to guarantee that the discrete reaction-diffusion system is bounded for all time. We consider reaction-diffusion systems on a 1D domain with homogeneous Neumann boundary conditions and nonnegative initial data and solutions. We define a Lyapunov-like function and show that its existence guarantees that the discrete reaction-diffusion system is bounded. These results are considered in the context of four example systems for which Lyapunov-like functions can and cannot be found.

WEAK SINDy: GALERKIN-BASED DATA-DRIVEN MODEL SELECTION.

We present a novel weak formulation and discretization for discovering governing equations from noisy measurement data. This method of learning differential equations from data fits into a new class of algorithms that replace pointwise derivative approximations with linear transformations and variance reduction techniques. Compared to the standard SINDy algorithm presented in [S. L. Brunton, J. L. Proctor, and J. N. Kutz, Proc. Natl. Acad. Sci. USA, 113 (2016), pp. 3932-3937], our so-called weak SINDy (WSINDy) algorithm allows for reliable model identification from data with large noise (often with ratios greater than 0.1) and reduces the error in the recovered coefficients to enable accurate prediction. Moreover, the coefficient error scales linearly with the noise level, leading to high-accuracy recovery in the low-noise regime. Altogether, WSINDy combines the simplicity and efficiency of the SINDy algorithm with the natural noise reduction of integration, as demonstrated in [H. Schaeffer and S. G. McCalla, Phys. Rev. E, 96 (2017), 023302], to arrive at a robust and accurate method of sparse recovery.

Competitive exclusion in a DAE model for microbial electrolysis cells

Microbial electrolysis cells (MECs) are devices that employ electroactive bacteria to perform extracellular electron transfer, enabling hydrogen generation from biodegradable substrates. In our previous work, we developed and analyzed a differential-algebraic equation (DAE) model for MECs. The model resembles a chemostat or continuous stirred tank reactor (CSTR). It consists of ordinary differential equations for concentrations of substrate, microorganisms, and an extracellular mediator involved in electron transfer. There is also an algebraic constraint for electric current and hydrogen production. Our goal is to determine the outcome of competition between methanogenic archaea and electroactive bacteria, because only the latter contribute to electric current and the resulting hydrogen production. We investigate asymptotic stability in two industrially relevant versions of the model. An important aspect of many chemostat models is the principle of competitive exclusion. This states that only microbes which grow at the lowest substrate concentration will survive as t → ∞.We show that if methanogens can grow at the lowest substrate concentration, then the equilibrium corresponding to competitive exclusion by methanogens is globally asymptotically stable. The analogous result for electroactive bacteria is not necessarily true. In fact we show that local asymptotic stability of competitive exclusion by electroactive bacteria is not guaranteed, even in a simplified version of the model. In this case, even if electroactive bacteria can grow at the lowest substrate concentration, a few additional conditions are required to guarantee local asymptotic stability. We provide numerical simulations supporting these arguments. Our results suggest operating conditions that are most conducive to success of electroactive bacteria and the resulting current and hydrogen production in MECs. This will help identify when producing methane or electricity and hydrogen is favored.

Life cycle of a cyanobacterial carboxysome.

Carboxysomes, prototypical bacterial microcompartments (BMCs) found in cyanobacteria, are large (~1 GDa) and essential protein complexes that enhance CO2 fixation. While carboxysome biogenesis has been elucidated, the activity dynamics, lifetime, and degradation of these structures have not been investigated, owing to the inability of tracking individual BMCs over time in vivo. We have developed a fluorescence-imaging platform to simultaneously measure carboxysome number, position, and activity over time in a growing cyanobacterial population, allowing individual carboxysomes to be clustered on the basis of activity and spatial dynamics. We have demonstrated both BMC degradation, characterized by abrupt activity loss followed by polar recruitment of the deactivated complex, and a subclass of ultraproductive carboxysomes. Together, our results reveal the BMC life cycle after biogenesis and describe the first method for measuring activity of single BMCs in vivo.

Mechanical regulation of photosynthesis in cyanobacteria.

Photosynthetic organisms regulate their responses to many diverse stimuli in an effort to balance light harvesting with utilizable light energy for carbon fixation and growth (source-sink regulation). This balance is critical to prevent the formation of reactive oxygen species that can lead to cell death. However, investigating the molecular mechanisms that underlie the regulation of photosynthesis in cyanobacteria using ensemble-based measurements remains a challenge due to population heterogeneity. Here, to address this problem, we used long-term quantitative time-lapse fluorescence microscopy, transmission electron microscopy, mathematical modelling and genetic manipulation to visualize and analyse the growth and subcellular dynamics of individual wild-type and mutant cyanobacterial cells over multiple generations. We reveal that mechanical confinement of actively growing Synechococcus sp. PCC 7002 cells leads to the physical disassociation of phycobilisomes and energetic decoupling from the photosynthetic reaction centres. We suggest that the mechanical regulation of photosynthesis is a critical failsafe that prevents cell expansion when light and nutrients are plentiful, but when space is limiting. These results imply that cyanobacteria must convert a fraction of the available light energy into mechanical energy to overcome frictional forces in the environment, providing insight into the regulation of photosynthesis and how microorganisms navigate their physical environment.

Medial septum activation produces opposite effects on dopamine neuron activity in the ventral tegmental area and substantia nigra in MAM vs. normal rats.

The medial septum (MS) differentially impacts midbrain dopamine (DA) neuron activity via the ventral hippocampus, a region implicated in DA-related disorders. However, whether MS regulation of ventral tegmental area (VTA) and substantia nigra pars compacta (SNc) is disrupted in a developmental disruption model of schizophrenia is unknown. Male Sprague-Dawley rats were exposed at gestational day 17 to methylazoxymethanol (MAM) or saline. As adults, NMDA (0.75 µg/0.2 µL) was infused into the MS, and either DA neuron activity in the VTA and SNc (7-9 anesthetized rats per group) or amphetamine-induced hyperlocomotion (AIH, 11-13 rats per group) was measured. MS activation produced a 58% increase in the number of spontaneously active DA neurons in VTA and a 37% decrease in SNc in saline rats. However, MS activation produced opposite effects on DA population activity in MAM rats, decreasing VTA DA activity by 51% and increasing SNc DA activity by 47%. MS activation also increased AIH by 113% in MAM rats, opposite of what is seen in intact rats. The effect in behavioral output may be due to disrupted GABAergic regulation of SNc as bicuculline infusion into vSub, which selectively prevented the MS activation-induced decrease in SNc DA activity in intact rats, prevented the increase in AIH and SNc DA activity in MAM rats. These findings demonstrate that the regulation of midbrain DA neurons by the MS is disrupted in this well-validated animal model, suggesting that it could be a potential locus for pharmacological intervention in disorders such as schizophrenia.

Medial septum differentially regulates dopamine neuron activity in the rat ventral tegmental area and substantia nigra via distinct pathways.

The medial septum (MS) impacts hippocampal activity and the hippocampus, in turn, regulates midbrain dopamine (DA) neuron activity. However, it remains to be determined how MS activation impacts midbrain DA activity. This question was addressed by infusing NMDA (0.75 µg/0.2 µL) into the medial septum of anesthetized male Sprague-Dawley rats and recording dopamine neuron activity in the ventral tegmental area (VTA) and substantia nigra pars compacta (SNc). MS activation increased (71%) the number of spontaneously active DA neurons in the VTA, and decreased (40%) the number of active DA neurons in the SNc. Effects in both the VTA and SNc required the ventral subiculum, but were differentially dependent on cholinergic and GABAergic mechanisms within the vSub and rostral and caudal subregions of the ventral pallidum, respectively. MS activation also decreased amphetamine-induced locomotor behavior, which was dependent on GABAergic inputs to the hippocampus. These findings demonstrate that the MS differentially regulates meso-striatal DA transmission via distinct pathways.

A numerical framework for computing steady states of structured population models and their stability.

Structured population models are a class of general evolution equations which are widely used in the study of biological systems. Many theoretical methods are available for establishing existence and stability of steady states of general evolution equations. However, except for very special cases, finding an analytical form of stationary solutions for evolution equations is a challenging task. In the present paper, we develop a numerical framework for computing approximations to stationary solutions of general evolution equations, which can also be used to produce approximate existence and stability regions for steady states. In particular, we use the Trotter-Kato Theorem to approximate the infinitesimal generator of an evolution equation on a finite dimensional space, which in turn reduces the evolution equation into a system of ordinary differential equations. Consequently, we approximate and study the asymptotic behavior of stationary solutions. We illustrate the convergence of our numerical framework by applying it to a linear Sinko-Streifer structured population model for which the exact form of the steady state is known. To further illustrate the utility of our approach, we apply our framework to nonlinear population balance equation, which is an extension of well-known Smoluchowski coagulation-fragmentation model to biological populations. We also demonstrate that our numerical framework can be used to gain insight about the theoretical stability of the stationary solutions of the evolution equations. Furthermore, the open source Python program that we have developed for our numerical simulations is freely available from our GitHub repository (github.com/MathBioCU).

A numerical framework for computing steady states of structured population models and their stability.

Structured population models are a class of general evolution equations which are widely used in the study of biological systems. Many theoretical methods are available for establishing existence and stability of steady states of general evolution equations. However, except for very special cases, finding an analytical form of stationary solutions for evolution equations is a challenging task. In the present paper, we develop a numerical framework for computing approximations to stationary solutions of general evolution equations, which can also be used to produce approximate existence and stability regions for steady states. In particular, we use the Trotter-Kato Theorem to approximate the infinitesimal generator of an evolution equation on a finite dimensional space, which in turn reduces the evolution equation into a system of ordinary differential equations. Consequently, we approximate and study the asymptotic behavior of stationary solutions. We illustrate the convergence of our numerical framework by applying it to a linear Sinko-Streifer structured population model for which the exact form of the steady state is known. To further illustrate the utility of our approach, we apply our framework to nonlinear population balance equation, which is an extension of well-known Smoluchowski coagulation-fragmentation model to biological populations. We also demonstrate that our numerical framework can be used to gain insight about the theoretical stability of the stationary solutions of the evolution equations. Furthermore, the open source Python program that we have developed for our numerical simulations is freely available from our GitHub repository (github.com/MathBioCU).