Comments on Mathematical Aspects of the Biró–Néda Model

Abstract

We address two mathematical aspects of the Biró–Néda dynamical model, recently applied in the statistical analysis of several and varied complex phenomena. First, we show that a given implicit assumption ceases to be valid outside the most simple and common cases, and we analyze the consequences thereof, in what the formulation of the model and probability conservation is concerned. Second, we revisit the transient behavior in the case of a constant reset rate and a constant or linear growth rate, improving on a previous analysis by including more general initial conditions.

1. Introduction

In recent years, Biró, Néda, and collaborators have successfully applied a simple dynamical model to address the behavior of complex systems of a varied nature, including high energy physics systems, ecology, population distribution, scientific citations, and social media popularity [1,2,3,4]. In addition, particularly promising are more recent applications to income and wealth distributions [5,6,7]. These latter applications contribute to an ongoing effort to model social and economic phenomena by means of simple mathematical models, often inspired in statistical physics and stochastic processes (see e.g., [8,9,10]). Going beyond the standard master equation for diffusion in a discrete system, the authors intend to explore unidirectional growth processes in order to model the dynamics of complex systems, including open ones. These processes explicitly break the detailed balance condition, and the unidirectional growth is supplemented with a mechanism designed to ensure the existence of interesting stationary distributions [1,2].

Specifically, the model concerns a probability distribution $\{P_n,n\ge 0\}$ evolving in time according to the following differential equations:

$$\dot{P_n}=-({\gamma }_n+{\mu }_n)P_n+{\mu }_{n-1}{P}_{n-1},n\ge 1,$$

(1)

$${\dot{P}}_0=\sum _{n=0}^{\infty }{\gamma }_n{P}_n-({\gamma }_0+{\mu }_0)P_0.$$

(2)

Within this class of models, each particular case is parametrized by the sequences $\left\{{\mu }_n\right\}$ and $\left\{{\gamma }_n\right\}$. One can recognize in (1) the differential equations for the probability of each state in a continuous-time Markov chain [11,12]. In particular, $\left\{{\mu }_n\right\}$ characterizes the dynamics of a pure birth process [11,12]. However, the one-step transitions are supplemented by nonlocal transitions to the state $n=0$, modifying the dynamics in order to allow for nontrivial stationary solutions.

In fact, the time-independent version of (1), i.e., with $\dot{P_n}=0$, for all $n\ge 1$, can be used by iteration as a generator for probability distributions, starting from some $Q_0>0$. Several interesting distributions have been considered and discussed in applications, for various $\left\{{\gamma }_n\right\}$ and $\left\{{\mu }_n\right\}$. In particular, for constant ${\gamma }_n$, the standard exponential distribution $Q_n=Q_0{e}^{-\beta n}$ emerges for the case of constant ${\mu }_n$, whereas a linear growth ${\mu }_n=\sigma (n+b)$ produces the so-called Waring distribution [13]. Other distributions of interest are generated by other sequences, with constant and nonconstant ${\gamma }_n$ and a variety of behaviors for ${\mu }_n$, including power-law and exponential growth (see Refs. [1,2], and also [14] for the emergence of distributions in a related context).

The purpose of the current article is twofold. In Section 4, we address the behavior of the transient, time-dependent solutions of the Biró–Néda model for the two most simple and common cases, namely with constant ${\gamma }_n$ and constant or linear ${\mu }_n$. In this respect, we develop the analysis first put forward in Ref. [15], hopefully contributing to the effort of exploring the behavior of these models, and in particular to improve the understanding of the approach to stationarity.

In the first part of the article, we point out and address technical challenges faced by the Biró–Néda model, outside the simple cases mentioned above. First, we will argue in Section 2 that the current formulation of the model requires modification, in order to accommodate the cases such that the sequence of values ${\mu }_n/{\gamma }_n$ grow faster than n. In fact, for those cases, we will show that the model does not actually admit (nontrivial) stationary solutions. The above-mentioned iteration process is of course present, however the only stationary solution compatible with (2) gives $Q_0=0$. On the technical side, we disprove a claim made in Ref. [1], showing that it fails for those fast sequences ${\mu }_n/{\gamma }_n$.

As a second contribution, we discuss the solution to the above issue. We will see in Section 3 that there is a natural and simple modification of Equation (2) compatible with the stationary solutions, however it comes at a cost. In fact, we will argue that probability is only conserved close to stationarity, which in this case means that it can only be conserved starting from initial distributions that are quite different from typical initial states. This consequence seems to be unavoidable, since the stationary solution introduces a nonzero boundary term in the equation for probability conservation.

It should again be said that no such issues affect the most simple and most common applications discussed by Biró, Néda, and collaborators, which most frequently involve simple situations for which the above-mentioned claim made in Ref. [1] is actually valid. Nevertheless, the departure from those simple cases is relevant, and our remarks shed light on the limits of applicability of the model, hopefully contributing to improve the mathematical formulation and the mathematical status of the model itself.

2. Stationary Solution

In this section we will disprove a general claim made in Ref. [1] and show that, for the cases where that claim is not valid, the current formulation of the Biró–Néda model [1,2] does not admit a nontrivial stationary solution.

The Biró–Néda model consists of the system of Equations (1) and (2) for the evolution of a probability distribution $\{P_n\left(t\right),n\ge 0\}$, where the dot stands for time derivative, and time t takes values in $[0,\infty [$. Following Refs. [1,2], we introduce parameters ${\lambda }_n$:

$${\lambda }_n={\mu }_n+{\gamma }_n,n\ge 0,$$

(3)

where all values ${\mu }_n$ and ${\gamma }_n$ are positive. The sequence of values ${\gamma }_n$ is referred to by the authors of Refs. [1,2] as the loss rate, or reset rate, and ${\mu }_n$ as the growth rate.

Let us then suppose that a stationary solution of the model exists, i.e., such that ${\dot{P}}_n=0$, for all t, for all $n\ge 0$. Let $P_n\left(t\right)=Q_n$ be that solution, where $\{Q_n,n\ge 0\}$ is some sequence of real numbers satisfying:

$$0=-{\lambda }_n{Q}_n+{\mu }_{n-1}{Q}_{n-1},\forall n\ge 1,$$

(4)

and a corresponding condition coming from Equation (2), which is:

$$\sum _{n=0}^{\infty }{\gamma }_n{Q}_n=Q_0{\lambda }_0.$$

(5)

On the other hand, one can easily obtain the expression for $Q_n$, for all $n\ge 1$, in terms of $Q_0$, by a standard iteration typical of these systems [12]. Starting from $Q_0$, the iteration of system (4) gives:

$$Q_n=Q_0\frac{{\mu }_0{\mu }_1\dots {\mu }_{n-1}}{{\lambda }_1{\lambda }_2\dots {\lambda }_n},n\ge 1.$$

(6)

Following Ref. [1], let us rewrite (6) in the form (valid also for $n=0$):

$$Q_n=\frac{Q_0{\lambda }_0}{{\gamma }_n}{r}_n\prod _{k=0}^n\frac{1}{1+r_k},\forall n\ge 0,\mathrm{where}{r}_k:=\frac{{\gamma }_k}{{\mu }_k}.$$

(7)

Again as in Ref. [1], we introduce also a special notation $S_0$ for the following infinite sum:

$$S_0=\sum _{n=0}^{\infty }{Z}_n,$$

(8)

where we have defined:

$$Z_n:=r_n\prod _{k=0}^n\frac{1}{1+r_k}.$$

(9)

Note that $S_0$ depends exclusively on the sequence $\{r_n={\gamma }_n/{\mu }_n\}$.

It follows immediately from (7) that:

$$\sum _{n=0}^{\infty }{\gamma }_n{Q}_n=Q_0{\lambda }_0{S}_0.$$

(10)

Thus, the authors of Refs. [1,2] have two distinct expressions for ${\sum }^{\infty }{\gamma }_n{Q}_n$: Equation (10) above and Equation (5), obtained from the evolution equation for $P_0$ (2). To reconcile these two expressions, it is claimed in Ref. [1] that $S_0$ always takes the value $S_0=1$, for all sequences $\left\{r_n\right\}$. We show next that this is not the case. In fact, and although $S_0$ is indeed equal to one for the most common and simple type of sequences $\left\{{\mu }_n\right\}$, $\left\{{\gamma }_n\right\}$ considered in Refs. [1,2], it turns out that $S_0$ is strictly smaller than one for more general sequences $\left\{r_n\right\}$.

Lemma 1.

$$\sum _{n=0}^N{Z}_n=1-Z_N/r_N.$$

(11)

Proof. 

Follows from the sum of Equation (4), or alternatively, directly from (9), by working out the expression for ${\sum }^N{Z}_n-1$. □

Proposition 1.

The infinite sum $S_0$ exists for all (positive) sequences $\left\{r_n\right\}$, with $S_0\le 1$. The equality $S_0=1$ is achieved if and only if the sum ${\sum }_{n=0}^{N}ln(1+r_n)$ diverges, when $N\to \infty$.

Proof. 

From Lemma 1, it follows that:

$$S_0=1-\underset{N\to \infty }{lim}\prod _{n=0}^N\frac{1}{1+r_n}.$$

(12)

The latter limit always exists, since the sequence is clearly decreasing and bounded from below. On the other hand, applying logarithm to the product, it follows that the limit of ${\prod }^N\frac{1}{1+r_n}$ is zero if and only if the sum ${\sum }^{N}ln(1+r_n)$ diverges. □

Note that $S_0=1$ is ensured for all sequences $\left\{r_n\right\}$ that do not tend to zero, and also for sequences $\left\{r_n\right\}$ that behave like $1/n$, for large n. However, $S_0$ fails to be equal to one for sequences $\left\{r_n\right\}$ that behave, for large n, as $r_n\sim n^{-1-ϵ}$, with $ϵ>0$, since the infinite sum ${\sum }^{\infty }ln(1+r_n)$ behaves in that case (apart from some finite term) as ${\sum }^{\infty }{n}^{-1-ϵ}$, which converges. Thus, $S_0=1$ is indeed satisfied for the simplest cases considered in Refs. [1,2,5,6,7,15], and in particular for constant sequences $\left\{{\gamma }_n\right\}$, with $\left\{{\mu }_n\right\}$ being also constant or linearly growing with n, but it fails even for constant $\left\{{\gamma }_n\right\}$ sequences if $\left\{{\mu }_n\right\}$ grows like $n^{1+ϵ}$.

If $S_0$ is not equal to one, then the Biró–Néda model simply does not admit a nontrivial stationary solution. In fact, the only solution to the joint Equations (5) and (10) is in that case $Q_0=0$, since ${\lambda }_0\ne 0$. This of course implies $Q_n=0$, for all n.

The source of the mismatch between expressions (5) and (10) present in Refs. [1,2] becomes clear when considering the derivative of the total probability. Let us then apply the time derivative to the infinite sum ${\sum }^{\infty }{P}_n$, taking into account (1) (and assuming that the derivative and limit of the sum commute). We obtain:

$$\frac{d}{dt}\sum _{n=0}^{\infty }{P}_n={\dot{P}}_0+\underset{N\to \infty }{lim}\sum _{n=1}^N{\dot{P}}_n={\dot{P}}_0-\underset{N\to \infty }{lim}\sum _{n=1}^N{\lambda }_n{P}_n+\underset{N\to \infty }{lim}\sum _{n=1}^N{\mu }_{n-1}{P}_{n-1}.$$

(13)

A rearrangement of the finite sums then gives:

$$\frac{d}{dt}\sum _{n=0}^{\infty }{P}_n={\dot{P}}_0-\underset{N\to \infty }{lim}\sum _{n=0}^N{\gamma }_n{P}_n+({\mu }_0+{\gamma }_0)P_0-\underset{N\to \infty }{lim}{\mu }_N{P}_N.$$

(14)

Typically, one would like to neglect the last term on the right-hand side, and the requirement of probability conservation would lead to Equation (2). However, this is incompatible with the possibility of a nontrivial stationary solution, if $S_0\ne 1$. In fact, for a distribution $\left\{Q_n\right\}$ generated from $Q_0$ by means of the system (4), it follows from (7) and (12) precisely that ${lim}_{N\to \infty }{\mu }_N{Q}_N={\lambda }_0{Q}_0(1-S_0)$, and thus ${lim}_{N\to \infty }{\mu }_N{Q}_N$ is zero for a nontrivial stationary solution if and only if $S_0=1$.

Equation (2) therefore does not admit (nontrivial) stationary solutions when $S_0\ne 1$. In order to include stationary solutions in the model, a modified version of Equation (2) is required, which we discuss in the next section.

To conclude the present section, let us remark that, regardless of the differential equation for $P_0$, the (nontrivial) stationary solution to be defined by (7) is not necessarily normalizable for arbitrary sequences $\left\{{\gamma }_n\right\}$, $\left\{{\mu }_n\right\}$. In fact, taking ${\gamma }_n={\mu }_n=e^{-n}$ leads to $Q_n\propto {(e/2)}^n$, which is clearly not summable. On the other hand, considering in particular situations with $S_0\ne 1$, normalizability is definitely ensured for the cases of fast growth rates described e.g., in Refs. [1,2]. In fact, it follows from $S_0\ne 1$ that $Q_n$ behaves asymptoticaly as $1/{\mu }_n$, and is therefore normalizable for all sequences ${\mu }_n$ growing faster than n

3. Dynamics and Probability Conservation

We will now obtain and discuss the generalization of Equation (2) to the cases where $S_0\ne 1$, allowing for nontrivial stationary solutions.

3.1. General Case

The simple way out of the conflict between Equations (5) and (10) is to include the boundary term ${lim}_{n\to \infty }{\mu }_n{Q}_n=Q_0{\lambda }_0(1-S_0)$ as an extra nonhomogeneous contribution in Equation (2), or alternatively to redefine the parameters ${\gamma }_n$ appearing in (2) (and only those, Equation (1) is left unchanged). Let us start with a free parameter $R>0$, define:

$${\tilde{\gamma }}_n={\gamma }_n+R{\lambda }_0(1-S_0),$$

(15)

and adopt the following equation for ${\dot{P}}_0$:

$${\dot{P}}_0=\sum _{n=0}^{\infty }{\tilde{\gamma }}_n{P}_n-{\lambda }_0{P}_0.$$

(16)

For stationary solutions we now obtain from (16):

$$\sum _{n=0}^{\infty }{\gamma }_n{Q}_n={\lambda }_0{Q}_0-R{\lambda }_0(1-S_0)\sum _{n=0}^{\infty }{Q}_n,$$

(17)

and a solution compatible with (10) emerges, with $Q_0=R$ and ${\sum }_{n=0}^{\infty }{Q}_n=1$. The value of R, and therefore of $Q_0$, is thus fixed by normalization.

There is however a price to be paid, in terms of probability conservation, for including stationary solutions in the model, with $S_0\ne 1$, precisely due to the nonzero boundary term. To see this, let us go back to the general Equation (14) and introduce Equation (16), to obtain:

$$\frac{d}{dt}\sum _{n=0}^{\infty }{P}_n=Q_0{\lambda }_0(1-S_0)\sum _{n=0}^{\infty }{P}_n-\underset{N\to \infty }{lim}{\mu }_N{P}_N.$$

(18)

(Alternatively, adding the boundary term simply as a nonhomogeneous contribution in (2) would have the only effect of removing the sum ${\sum }^{\infty }{P}_n$ from the right-hand side).

Equation (18) indicates that probability can only be conserved for initial probability distributions $\left\{P_n\left(0\right)\right\}$ such that ${lim}_{n\to \infty }{\mu }_n{P}_n\left(0\right)=Q_0{\lambda }_0(1-S_0)$. When the boundary term contribution is absent, the typical initial distributions, such that $P_n\left(0\right)=0$ for all n greater than a certain $n_0$, would ensure a null right-hand side in (18). This is no longer the case if $S_0\ne 1$. According to (18), probability is only conserved close to stationarity, e.g., for distributions with the same asymptotic behavior as the stationary solution. The usefulness of this model therefore seems to be restricted in cases where $S_0\ne 1$, since the states with $P_n\left(0\right)=0$ for all $n>n_0$ are typical initial states.

3.2. Constant Reset Rate

As preparation for the next section, and also for its own intrinsic interest, we now consider the special case of constant $\left\{{\gamma }_n\right\}$ sequences, i.e., with ${\gamma }_n=\gamma$, for all n.

Considering Equation (16), it is certainly tempting to simply impose ${\sum }^{\infty }{P}_n=1$ for all t and leave the remaining unchanged. However, since probability conservation is itself an issue (for $S_0\ne 1$), it is safer to start afresh with a nonhomogeneous equation for ${\dot{P}}_0$ of the form:

$${\dot{P}}_0=\gamma -{\lambda }_0{P}_0+K,$$

(19)

where K is some constant, and see where it leads. Equation (19) now gives $Q_0$ in terms of K, and compatibility with (10) and normalization of $\left\{Q_n\right\}$ fixes K as $K=\gamma /S_0-\gamma$. We therefore obtain:

$${\dot{P}}_0=-{\lambda }_0{P}_0+\frac{\gamma }{S_0},$$

(20)

whereas (18) is replaced with:

$$\frac{d}{dt}\sum _{n=0}^{\infty }{P}_n=\gamma \left(1-\sum _{n=0}^{\infty }{P}_n\right)+Q_0{\lambda }_0(1-S_0)-\underset{N\to \infty }{lim}{\mu }_N{P}_N.$$

(21)

Thus, with ${\sum }^{\infty }{P}_n=1$, both Equations (20) and (21) coincide with those obtained from (16)–(18).

It is worthwhile mentioning the following. As we have seen, for $S_0\ne 1$, Equation (2) is fully incompatible with (nontrivial) stationary solutions. This is not the case with its reduced counterpart, e.g., an equation of the form ${\dot{P}}_0=-{\lambda }_0{P}_0+\gamma$, since (5) would be replaced by $\gamma =Q_0{\lambda }_0$, which is of a different nature, and not fully incompatible with (10). In fact, a stationary solution exists, with ${\sum }_{n=0}^{\infty }{Q}_n=S_0$. However, this does not mean that an unnormalized stationary solution for Equation (2) exists, as it does not. Note also that the reduction from Equation (2) to its reduced counterpart would be inconsistent if $S_0\ne 1$, since it would be done assuming ${\sum }_{n=0}^{\infty }{P}_n=1$, whereas the reduced equation gives ${\sum }_{n=0}^{\infty }{Q}_n=S_0$.

Returning to Equation (21), the situation regarding probability conservation remains as discussed in Section 3.1 above. Examples leading to $S_0\ne 1$ are ${\mu }_n=n^s+1$, with $s>1$. In this case, any initial probability distribution decaying faster than $1/n^s$ fails to produce a null right-hand side in (21). For bounded or linearly growing sequences $\left\{{\mu }_n\right\}$, on the other hand, the model formed by system (1), with ${\gamma }_n=\gamma$, and Equation (20) (with $S_0=1$) is fully consistent to describe the time evolution of a probability distribution $\left\{P_n\right\}$.

4. Behavior of the Transient Solution

In this section, we depart from foundational and formulation issues, and address the behavior of the time-dependent solutions of the Biró–Néda model for two special cases. In particular, we consider sequences $\left\{{\gamma }_n\right\}$ that are constant, and sequences $\left\{{\mu }_n\right\}$ that are either constant or growing linearly with n. In both cases the value of $S_0$ (8) is equal to one, no boundary issues such as those discussed in Section 2 and Section 3 exist, and the model is fully consistent, as we have just seen.

This issue was originally addressed in Ref. [15]. In that work, the general solutions of the differential equations were written down, for the two considered cases, and an attempt was made to study the behavior of the transient, i.e., time-dependent, solution. Particular attention was paid to the monotonicity, or lack thereof, of the time-dependent solution when approaching the stationary solution. Although the general solutions presented in Ref. [15] are absolutely correct, the analysis of their behavior was restricted to a very special case, embodied by very particular initial conditions. We show here that the simple behavior of the solutions discussed in Ref. [15] is far from generic.

Let us then fix ${\gamma }_n=\gamma$, for all $n\ge 0$, and let $\left\{P_n\right\}$ be a sequence of variables satisfying Equations (1) and (20) with $S_0=1$, for sequences $\left\{{\mu }_n\right\}$ that grow, at most, linearly with n.

Let us write, for arbitrary n:

$$P_n\left(t\right)={\Delta }_n\left(t\right)+Q_n.$$

(22)

Since, for each fixed n, the set $\{Q_0,\cdots ,Q_n\}$ is a solution of the nonhomogeneous system of differential equations determined by (1) and (20) for the variables $\{P_0,\cdots ,P_n\}$, it follows from general arguments that the variables ${\Delta }_n$ are solutions of the corresponding homogeneous system, i.e.,

$$\begin{array}{c}\hfill {\dot{\Delta }}_0=-{\lambda }_0{\Delta }_0,\\ \hfill {\dot{\Delta }}_n=-{\lambda }_n{\Delta }_n+{\mu }_{n-1}{\Delta }_{n-1},n\ge 1.\end{array}$$

(23)

Given that all values ${\lambda }_n$, $n\ge 0$, are positive, the functions ${\Delta }_n$ tend to zero as time tends to infinity, with $P_n$ approaching its stationary value $Q_n$. The detailed behavior of the functions ${\Delta }_n$ was analyzed in Ref. [15] only for the particular initial conditions $P_0\left(0\right)=1$, $P_n\left(0\right)=0$, for all $n>0$. The authors observed only one stationary point for each function ${\Delta }_n$ and concluded (analytically in the case of constant ${\mu }_n$ and based on numerical evidence in the other case) that the functions ${\Delta }_n$ are monotonic after a critical time $t_c$ that depends on n. We show next that the behavior of the functions ${\Delta }_n$ is much more involved, in general. In fact, depending on the initial conditions, each function ${\Delta }_n$ can have up to n stationary points, with the position of the last one being arbitrary. Thus, while it is true that all functions ${\Delta }_n$ are eventually monotonic, the time after which monotonicity is achieved can be arbitrarily large, for any given n.

4.1. Constant Growth Rate

Let us then fix ${\mu }_n=\mu$ for all $n\ge 0$, besides ${\gamma }_n=\gamma$. As shown in Ref. [15], the solutions ${\Delta }_n$ to the homogeneous system determined by (23) can be written in terms of the initial values ${\Delta }_j\left(0\right)$, $j=0,\cdots ,n$ in the form:

$${\Delta }_n\left(t\right)=e^{-\lambda t}{f}_n\left(t\right),\lambda =\gamma +\mu ,$$

(24)

where $f_n$ is a degree-n polynomial:

$$f_n\left(t\right)=\sum _{k=0}^n{\Delta }_{n-k}\left(0\right)\frac{{\left(\mu t\right)}^k}{k!}.$$

(25)

The derivative of ${\Delta }_n$ is easily obtained:

$${\dot{\Delta }}_n=e^{-\lambda t}({\dot{f}}_n-\lambda f_n).$$

(26)

Thus, the stationary points of ${\Delta }_n$ coincide with the positive roots of the function ${\dot{f}}_n-\lambda f_n$, which is again a polynomial of degree n. In particular:

$${\dot{f}}_n-\lambda f_n=\sum _{k=0}^{n-1}\left[\mu {\Delta }_{n-k-1}\left(0\right)-\lambda {\Delta }_{n-k}\left(0\right)\right]\frac{x^k}{k!}-\lambda {\Delta }_0\left(0\right)\frac{x^n}{n!},x=\mu t.$$

(27)

Thus, each ${\Delta }_n$ has at most n stationary points, and is therefore monotonic after some point. However, for any given pair $\mu$, $\lambda$, one can find initial conditions such that, for any given n, the polynomial ${\dot{f}}_n-\lambda f_n$ is fairly arbitrary, concerning both the number of positive roots and the position of its largest root. In fact, the following applies:

Proposition 2.

Let μ, λ, and n be given, and let $t_n$ denote the largest (positive) root of ${\dot{f}}_n-\lambda f_n$. Then, for any $M>1$, one can find initial values ${\Delta }_k\left(0\right)$, $\forall k\ge 0$, such that ${\dot{f}}_n-\lambda f_n$ has n positive roots and $t_n>M$.

Proof. 

Let $p_n\left(x\right)={\sum }_{k=0}^n{a}_k{x}^n$ denote an arbitrary degree-n polynomial. The equation ${\dot{f}}_n-\lambda f_n=p_n$ gives us a system of linear equations for the set ${\Delta }_0\left(0\right),\cdots ,{\Delta }_n\left(0\right)$, in terms of the set $a_0,\cdots ,a_n$, which can always be solved. Now choose $p_n$ such that it possesses n positive roots, which is certainly possible (choose any polynomial with n real roots and apply a translation). If the largest root of $p_n$ is not greater than the prescribed value $\mu M$, replace $p_n\left(x\right)$ with $p_n(x/\mathrm{\Lambda })$, for appropriate $\mathrm{\Lambda }>1$. The choice ${\dot{f}}_n-\lambda f_n=p_n$ thus gives us n positive roots and an arbitrarily large value of $t_n$. To prove the proposition, it remains to show that the initial values ${\Delta }_k\left(0\right)$, $k\le n$, can be chosen in a way that is compatible with a probability distribution, namely defined by $P_k={\Delta }_k+Q_k$. It is thus necessary that:

$$-Q_k\le {\Delta }_k\left(0\right)\le 1-Q_k,$$

(28)

and

$$\sum _{k=0}^n\left({\Delta }_k\left(0\right)+Q_n\right)\le 1,$$

(29)

which follow from $0\le P_k\le 1$, $\forall k$, and ${\sum }_{k=0}^n{P}_k\le 1$. This can be simply achieved (if necessary) by replacing the latter polynomial with $ϵp_n\left(x\right)$, for sufficiently small $ϵ>0$. The new polynomial has the same set of roots, and it is certainly possible to satisfy all the conditions (28) concerning $k\le n$. Note again that the relation between the $a_k$’s and the ${\Delta }_k\left(0\right)$’s is linear, and therefore the scaling $p_n\to ϵp_n\left(x\right)$ results in a scaling ${\Delta }_k\left(0\right)\to ϵ{\Delta }_k\left(0\right)$. All conditions (28) for $k\le n$ can be satisfied for an appropriate choice of $ϵ$, since the number of conditions is finite. As for condition (29), it can only be violated if ${\sum }_{k=0}^n{\Delta }_k\left(0\right)>0$, since ${\sum }_{k=0}^n{Q}_k<1$ is satisfied by construction. In that case, one can just further reduce the value of $ϵ$ above, until (29) is satisfied. Clearly, this last action preserves conditions (28), therefore we are done for $k\le n$. For the remaining initial values one can simply set $P_{n+1}\left(0\right)=1-{\sum }_{k=0}^n\left({\Delta }_k\left(0\right)+Q_n\right)$ and $P_k=0$, $\forall k>n+1$. □

4.2. Linear Growth Rate

When all values ${\lambda }_n$ in the sequence $\{{\lambda }_n={\gamma }_n+{\mu }_n\}$ are different, the general solution for the functions ${\Delta }_n$ takes the form [15]:

$${\Delta }_n\left(t\right)=\sum _{k=0}^n{C}_k{\alpha }_k^n{e}^{-{\lambda }_{k}t},$$

(30)

with

$${\alpha }_k^n=\frac{{\prod }_{m=k}^{n-1}{\mu }_m}{{\prod }_{m=k+1}^n({\lambda }_m-{\lambda }_k)}.$$

(31)

The coefficients $C_k$ and the initial conditions ${\Delta }_k\left(0\right)$ are related by:

$${\Delta }_n\left(0\right)=\sum _{k=0}^n{C}_k{\alpha }_k^n,$$

(32)

or by the inverse relation:

$$C_k=\sum _{j=0}^k\prod _{i=j}^{k-1}\frac{{\mu }_i}{{\lambda }_i-{\lambda }_k}{\Delta }_j\left(0\right).$$

(33)

The above expressions take a simpler form in the particular case considered in Ref. [15]. Let then:

$${\lambda }_n=\gamma +\sigma (n+1),\sigma >0,\forall n\ge 0.$$

(34)

It follows from (31) that the coefficients ${\alpha }_k^n$ are just the binomial coefficients, i.e.,

$${\alpha }_k^n=\left(\begin{array}{c}n\\ k\end{array}\right).$$

(35)

In addition, expression (30) simplifies to:

$${\Delta }_n\left(t\right)=e^{-(\gamma +\sigma )t}\sum _{k=0}^n{C}_k{\alpha }_k^n{\left(e^{-\sigma t}\right)}^k,$$

(36)

with ${\alpha }_k^n$ given by (35). Relation (33) in turn becomes:

$$C_k=\sum _{j=0}^k{(-1)}^{k-j}\left(\begin{array}{c}k\\ j\end{array}\right){\Delta }_j\left(0\right).$$

(37)

Therefore, again ${\Delta }_n\left(t\right)$ is of the form:

$${\Delta }_n=e^{-\lambda t}{g}_n,$$

(38)

with $\lambda =\mu +\gamma$, where $g_n$ is a polynomial, this time in the variable $y=e^{-\sigma t}$, i.e.,

$$g_n=\sum _{k=0}^n{C}_k{\alpha }_k^n{y}^k,y=e^{-\sigma t}.$$

(39)

The derivative is:

$${\dot{\Delta }}_n=e^{-\lambda t}({\dot{g}}_n-\lambda g_n),$$

(40)

and again the stationary points of ${\Delta }_n$ correspond to the positive roots of the polynomial ${\dot{g}}_n-\lambda g_n$, which now takes the form:

$${\dot{g}}_n-\lambda g_n=-\sum _{k=0}^n{C}_k{\alpha }_k^n(\lambda +k\sigma )y^k.$$

(41)

Obviously, the stationary points t of ${\Delta }_n$ are related to the roots y of ${\dot{g}}_n-\lambda g_n$ by $t=-\frac{1}{\sigma }ln\left(y\right)$.

The following analogue of proposition (2) can be proven, using essentially the same arguments, with the appropriate adaptations.

Proposition 3.

Let γ, σ, and n be given, and let $t_n$ denote the largest (positive) stationary point of ${\Delta }_n$ (36). Then, for any $M>1$, one can find initial values ${\Delta }_k\left(0\right)$, $\forall k\ge 0$, such that ${\Delta }_n$ has n positive stationary points and $t_n>M$.

Proof. 

Let $p_n\left(y\right)={\sum }_{k=0}^n{a}_k{y}^n$ denote an arbitrary degree-n polynomial. Again, the equation ${\dot{g}}_n-\lambda g_n=p_n$ gives us a system of linear equations for the set ${\Delta }_0\left(0\right),\cdots ,{\Delta }_n\left(0\right)$, in terms of the set $a_0,\cdots ,a_n$, which can always be solved, taking into account (32). Now choose $p_n$ such that it possesses n roots in the interval $]0,1]$. This is certainly possible, since one can pick any polynomial $p_n\left(y\right)$ with n positive roots and replace it by $p_n(y/\mathrm{\Lambda })$, for appropriate $\mathrm{\Lambda }<1$. Let $y_0$ be the smallest root of $p_n$. If $t_n=-\frac{1}{\sigma }ln\left(y_0\right)$ is not greater than the prescribed value M, then simply adjust the value of $\mathrm{\Lambda }$ above. The choice ${\dot{g}}_n-\lambda g_n=p_n$ thus gives us n-positive stationary points and an arbitrarily large value of $t_n$. To conclude the proof, it remains to show that the initial values ${\Delta }_k\left(0\right)$ can be chosen such that $P_k={\Delta }_k\left(0\right)+Q_k$ defines a normalized probability distribution, which can be done by exactly the same arguments as in proposition (2). □

5. Conclusions and Discussion

We have analyzed the Biró–Néda model [1,2], parametrized by the reset rate ${\gamma }_n$ and growth rate ${\mu }_n$. We have shown that, contrary to expectations, the model does not actually admit (nontrivial) stationary solutions, if the sequence of values ${\mu }_n/{\gamma }_n$ grows faster than n. This follows from the fact that the sum $S_0$ discussed in Ref. [1] fails to be equal to one in such cases. Interesting time-independent distributions are nevertheless associated with the sequences ${\mu }_n$ and ${\gamma }_n$, however they simply do not appear as solutions of the model. This is related to the fact that those distributions $Q_n$ are such that ${lim}_{n\to \infty }{\mu }_n{Q}_n$ is nonzero, introducing an unbalanced boundary term at infinity.

Keeping the spirit of the Biró–Néda model, and trying in particular to remain within the setting of unidirectional processes, we have proposed a modified version of the dynamical equation for ${\dot{P}}_0$ (2), allowing for those distributions to emerge as true stationary solutions. This effort is partly successful, in the sense that the solutions are included in the model. However, probability is only conserved close to stationarity. This consequence looks unavoidable, since the lack of backwards interactions seems to leave no other option but to include the solutions essentially by introducing a nonhomogeneous term in the evolution of $P_0$, which is consequently reflected as a boundary term in the equation for the conservation of probability. This is to be contrasted e.g., with birth and death processes, where exactly the same stationary solutions can be obtained, with no boundary issues and with unrestricted probability conservation.

It is perhaps not surprising that technical challenges are arising in the Biró–Néda model for fast growth. In fact, at least for constant ${\gamma }_n$, the large n behavior is expected to be dominated by the ${\mu }_n$ terms in Equation (1). This is essentially a pure birth process, and those are known to exhibit exotic behavior for a fast growth rate, e.g., explosions (see e.g., [16,17]) and probability loss [18], appearing under essentially the same mathematical conditions found in Proposition 1.

Concerning other technical points, one needs to ensure that $\left\{Q_n\right\}$ is normalizable (which we have also shown is not a priori guaranteed for all sequences ${\gamma }_n$ and ${\mu }_n$) and that the infinite sum ${\sum }^{\infty }{\gamma }_n{P}_n$ appearing in the differential equation for $P_0$ makes sense. In fact, the coupling of $P_0$ with an infinite number of variables $P_n$ seems technically challenging, and a fast decaying sequence ${\gamma }_n$ would be preferable for a more amenable formulation. This however conflicts with the normalizability of $\left\{Q_n\right\}$. Requiring ${\gamma }_n\in ]a,b[$ for all n, for some strictly positive a and $a<b<\infty$, might be sufficient to address both questions, however a more detailed study is probably required.

Finally, in Section 4, we focused on the time-dependent behavior of the functions ${\Delta }_n=P_n-Q_n$, in the amenable cases of a constant reset rate and constant or linear growth rate. Our analysis considerably expands the discussion of the transient behavior presented in Section 6 of Ref. [15]. The behavior seen in Ref. [15] is not generic, and in fact comes from the very particular initial conditions $P_n\left(0\right)=0$, for all $n>0$, which, taking into account Equation (1), effectively compel all but one of the stationary points of $P_n$ to be located at $t=0$. In general, we were able to show that ${\Delta }_n$, and therefore $P_n$, can possess n stationary points, the last one occurring at an arbitrarily long time. Thus, monotonic convergence can occur, for general initial conditions, only after an arbitrarily long time.