Passivity analysis of discrete-time genetic regulatory networks with reaction-diffusion coupling and delay-dependent stability criteria

1.   Introduction

Genetic regulatory networks, often abbreviated as GRNs, are complex systems that include several components that work together to regulate gene expression. These networks encompass a wide range of interactions, including those between DeoxyriboNucleic Acid (DNA) and proteins (such as transcription factors), RiboNucleic Acid (RNA) and proteins (like ribosomes and RNA-binding proteins), and even interactions among different types of RNA. GRNs play a pivotal role in understanding various biological processes and have garnered significant attention from researchers ^[1]. Over the last couple of years, advances in molecular-level biology research have highlighted the significance of genetic networks, which are inherently complex and need simplification. GRNs have successfully described biochemical reactions such as gene transcription, translation, and protein diffusion, similar to the achievements in artificial neural network modeling, which have resulted in a plethora of models. Various models have emerged, including the Bayesian model ^[2], the Boolean model ^[3,4], and the differential equation model ^[5,6]. In particular, the differential equation model is good at describing the continuous evolution of a system state over time, and usually can reflect the real dynamics of the system more accurately. However, continuous systems are not suitable for computer processing, so many studies have converted continuous systems into discrete systems ^[7,8]. In this context, incorporating delayed discrete-time GRNs becomes crucial for simulation and implementation ^[9]. Numerous excellent results in ^[10,11] have analyzed exponential stability, made synchronization control, designed guaranteed cost control, and updated estimation methods for discrete-time GRNs. Recently, some notable contributions have further enriched this field. By resorting to the sliding-mode control, Li et al. ^[12] effectually extended and modified the existing works on the stabilization problem for delayed uncertain semi-Markovian switching complex-valued networks. They proposed sufficient criteria for stochastic stability based on the generalized Dynkin's formula and Lyapunov stability theory. In addition to these works, Narayanan et al. ^[13] addressed the problem of Mittag-Leffler synchronization of T-S fuzzy fractional-order discrete-time complex-valued molecular models of mRNA and protein in regulatory mechanisms with two kinds of regulation functions. Wei et al. ^[14] addressed the dissipation synchronization issue of semi-Markovian jumping delayed neural networks under random deception attacks and explained the specific influence of deception attacks on event trigger time for the first time. Pandiselvi et al. ^[15] studied the approximation concern for the discrete-time stochastic GRNs with the leakage delays, distributed delays, and probabilistic measurement delays into the problem and modeled the robust H_∞ state estimator for a class of discrete-time stochastic GRNs. Zhao and Wu ^[16] contributed to the field by establishing new fixed/prescribed-time stability criteria for stochastic systems with time delay and designed delay-independent control mechanisms that are applicable to systems with unknown delays. These studies highlight the importance of addressing key issues in GRN control and synchronization. The study of discrete-time GRN with spatial coupling has gained significant attention due to its importance in understanding the spatial and temporal dynamics of gene expression. In recent years, discrete-time GRNs with spatial coupling have emerged as a crucial area of research in systems biology as they integrate temporal and spatial dynamics to provide a more comprehensive understanding of gene regulation processes. Unlike existing literature that limits its focus to discrete-time networks, Zhang and Li ^[17] explored the double effects of both discrete time and discrete spatial diffusions in switching complex dynamical networks. Their work is ground breaking in its consideration of discrete spatial diffusions, offering a robust theoretical and practical foundation for future research in this area. Liang and Wang ^[18] studied the passivity for coupled reaction-diffusion neural networks with multiple state couplings or spatial diffusion couplings by employing the proportional-derivative control method. This research provides valuable insights into the control and synchronization of networks with reaction-diffusion terms. Discrete-time GRNs with spatial coupling have also been applied to study gene regulatory networks in disease and aging. Avila et al. ^[19] delineated different types of GRNs and provided their biological interpretation, highlighting the potential of discrete-time GRNs with spatial coupling to uncover the underlying mechanisms of complex diseases and aging processes.

Passivity, a mathematical principle describing a system or function's ability to absorb or retain energy under specific conditions, is often linked to energy retention and dissipation in control theory. At its core, passivity delineates a system's capacity to absorb, dissipate, or retain energy without amplifying it under specified conditions, thereby ensuring stability and boundedness in its state space. A dynamic system is deemed passive if its energy does not increase within its state space ^[20]. Recently, the concept of dissipation has found widespread application in various fields, including chemistry, biology, complex networks, signal processing, chaos control, and social sciences. Consequently, exploring the passivity of gene regulatory networks is crucial due to its relevance in numerous research and analytical spheres.

In genetic regulation, biochemical reactions involve the diffusion of regulatory proteins and metabolites between the cytosol and nucleus. It is more accurate to consider reaction-diffusion models for genetic regulatory networks than to assume spatial homogeneity. In fact, it has been suggested that reaction-diffusion systems play an important role in the dynamics of GRNs ^[21,22]. Mathematical models without considering reaction-diffusion effects may lead to incorrect predictions of protein and mRNA concentrations. Therefore, it is essential to consider reaction-diffusion systems in the models of GRNs. However, there are few relevant articles, especially in the study of discrete-time GRNs with reaction-diffusion in plants. In discrete-time gene regulatory networks, the nature of reaction-diffusion reflects the coupling between adjacent spatial locations. Specifically, it represents the interaction between neighboring spatial locations due to their close proximity, which leads to mutual influence on each other's function. Based on this observation, the present study establishes a gene regulatory network model that incorporates the mutual coupling between spatial locations. The GRNs constructed through differential equations in ^[5,6] are continuous models that solely reflect the temporal changes in protein and mRNA concentrations. In contrast to the continuous models proposed in ^[5,6], the proposed discretized GRN model offers additional insights. It can not only capture temporal variations in protein and mRNA concentrations but also reveal their spatial variations through the inclusion of reaction-diffusion terms. Computer simulations demonstrate the effectiveness of the theoretical results and the importance of incorporating reaction-diffusion terms in GRN models.

The main difficulties and challenges of this paper stem from the fact that the diffusion terms in GRNs first appeared in models described by differential equations, with limited research on discrete-time GRNs. For electronic computers, continuous models are not directly describable. Therefore, the modeling of continuous GRNs becomes highly necessary. However, since spatial diffusion is characterized by the spatial divergence of mRNA and protein concentrations, defining the spatial diffusion terms of mRNA and proteins in discrete-time GRN models and proving that such models exhibit passivity is quite challenging. In this paper, we address these challenges by discretizing spatial divergence and leveraging Lyapunov functions and linear matrix inequalities to provide solutions.

The main contributions of this article are as follows:

1) The introduction of the concept of spatial diffusion into a discrete GRN, constructing a discretely coupled gene regulatory network with reaction diffusion.

2) By considering Brownian motion and time delay, the article provides a theorem for determining the passivity in a discretely GRN with Dirichlet boundary conditions and reaction diffusion.

For readability, the meanings of commonly used mathematical symbols involved in this paper are as follows:

$R^{N}$ means the $N$ dimensional real vector, and $R^{N \times N}$ stands for the sets of $N \times N$ real matrices. $\frac{\partial f(x)}{\partial x}$ denotes the partial derivative of $f(x)$ with respect to $x$. $\sum\limits_{i = 1}^{N} \gamma_{i} = \gamma_{1} + \gamma_{2} + \cdots + \gamma_{N}$, $\bigcup\limits_{i = 1}^{N} \alpha_{i} = \{ \alpha_{1}, \alpha_{2}, \alpha_{N}\}$. The superscript $\rm T$ stands for matrix transposition.

2.   Problem formulation

A typical nonlinear delayed genetic regulatory network based on a differential equation model with reaction-diffusion terms is described as ^[21]:

where

and $B = (b_{nn'})\in R^{N \times N}$, $b_{nn'}$ represents the interaction between transcription factors and genes. Specifically, $b_{nn'}$ can take one of three values: $\alpha_{nn'}$, 0, or $-\alpha_{nn'}$. Here, $\alpha_{nn'}$ represents the magnitude of the interaction occurring between transcription factor $n'$ and gene $n$. In detail, a positive $\alpha _{nn'}$ signifies that there is an activating influence, a zero value indicates the absence of any direct interaction, and a negative $\alpha_{nn'}$ signifies that there is a repressing influence.

In Eq (2.1), the matrix $B$ represents the interactions between genes in the GRN. The notations $m'_{n}(t, x)$ and $p'_{n}(t, x)$ $(n = 1, 2, \cdots, N)$ represent the concentrations of mRNA and protein at time $t$ for the $n$th node, respectively. The vector $x = (x_{1}, x_{2}, \cdots, x_{l})^{\rm T}$ belongs to the set $Q$. Specifically, $Q$ is defined as $\{x| |x_{\alpha}| \leq L_{\alpha}, \alpha = 1, 2, \cdots, l\}$ and represents a compact set in the real vector space $R^{l}$ with a smooth boundary $\partial Q$. The terms $u(t, x)$ and $v(t, x)$ stand for the disturbance input signals. The parameters $A$ and $C$ are the degradation rates of the mRNA and protein, respectively, while $D$ signifies the translation rate. The matrices $M_{\alpha } > 0$ and $M_{\alpha}^{*} > 0$ demonstrate the diffusion coefficient matrices. The regulatory function $f_{n}(x)$ is in the Hill form and captures the feedback regulation of protein on transcription; its mathematical form is described as

where $H_{n}$ is the Hill coefficient and $v_{n}$ is a positive constant.

Based on Eq (2.1), we construct a corresponding discrete gene regulatory network model with reaction diffusion terms, as shown in Figure 1.

where

$m_{n}(k, i)$ and $p_{n}(k, i)$ represent the concentration of mRNA and protein of the $n$th node at the time $k \in N$ and space $i$ respectively, where $i = 1, 2, \cdots, L$. $\Delta_{1}$ and $\Delta_{2}$ represent the reaction diffusion terms of concentrations of mRNA and protein, respectively. The discretization of the spatial diffusion term is obtained by discretizing the spatial diffusion term of the continuous gene regulatory network as shown in Eq (2.3), with the specific correspondence shown in the following equations:

In the modeling of gene regulatory networks, it is very important to choose a suitable discretization scheme. Different discretization schemes may have advantages in different application scenarios. For example, the forward difference and backward difference methods may be easier to implement or more efficient in some cases, although their accuracy is low. In our model, the discretization formula for the second-order partial derivative $\frac{\partial^{2} m'(t, x)}{\partial x^{2}_{k}}$ typically employs the central difference method for approximation because this method estimates the derivative of a point by calculating the values on both sides of the functio. This can provide a more accurate derivative approximation and is very important for understanding and predicting the behavior of gene regulatory networks. In addition, the central difference method can better maintain the physical meaning and mathematical properties of the system when dealing with the spatial diffusion term.

For the function $m'(t, x)$, the discretization formula for its second-order partial derivative is as follows:

where $h_{k}$ is the step size in the direction of $x_{k}$.

After the model is discretized, the spatial diffusion term also needs to be discretized. The term $m'(t, x)$ is transformed into $m(t, i)$, where $t$ corresponds to time, and $i$ corresponds to the spatial position $x$. The discretization process is carried out to facilitate computer simulation and computation. Clearly, in discrete systems, the step size is typically taken as 1, that is, $h_{k} = 1$, and Eq (2.5) can be simplified as follows:

If $M_{k}$ is set to $\theta$, then the reaction-diffusion term of mRNA concentration

in the discrete form can be obtained. Similarly, the reaction-diffusion term $\Delta_{2}$ for protein concentration can be obtained in the discrete form, and thus we arrive at Eq (2.8).

where $\theta$ and $\mu$ are positive coupling coefficients.

Once Eq (2.8) has been substituted into Eq (2.3), and the subsequent simplification process has been carried out, the final expression is obtained:

where $A = A'-2\theta$ and $C = C'-2\mu$.

Remark 2.1. To provide a concrete application example, we refer to the research by Zhang et al. ^[23], which effectively demonstrates the importance of combining regulatory networks with reaction-diffusion terms. Zhang et al. studied the oscillatory expression of Escherichia coli mediated by microRNA, incorporating both time delays and reaction-diffusion components. They demonstrated that the oscillatory expression of Escherichia coli is not only critically dependent on transcriptional and translational delays but also significantly influenced by the diffusion coefficients. This conclusion has been practically verified by numerous biological experiments and observations. They also found that if the diffusion coefficients of microRNA, mRNA, and protein are sufficiently small, non-uniform periodic oscillations can be predicted to occur, unless only spatially uniform periodic oscillations are exhibited. Thus, the application example provided by Zhang et al. ^[23] serves as a practical illustration of how the combination of regulatory networks and reaction-diffusion terms, as considered in our paper, can be applied to model and understand complex biological phenomena such as the oscillatory expression of Escherichia coli.

The current research focuses on the Dirichlet boundary condition, which is characterized by the following mathematical expressions:

Here, the symbol $\partial \Omega$ denotes the boundary of the domain under consideration.

Dirichlet boundary conditions provide a clear and physically meaningful way to define the behavior of the system at the boundaries, which is very important for drawing meaningful conclusions about the passivity and stability of the network. The use of Dirichlet boundary conditions in our model enables us to focus on the internal dynamics of the system instead of being dispersed by the complexities of boundary behavior. This simplification enables us to deduce the analysis results, which provide a deep understanding of the behavior of gene regulatory networks with spatial diffusion. The unique contribution of our research lies in applying these boundary conditions to solve the specific problem of analyzing the passivity of discrete gene regulatory networks with spatial diffusion.

Considering that the concentrations of mRNA and proteins are influenced by molecular Brownian motion, Eq (2.9) is rewritten as

where $\omega(k, i)$ represents a vector-form scalar Brownian motion with the following properties:

and $\sigma\Big(p(k, i), p\big(k-\tau(k), i\big)\Big)$ is the noise intensity matrix satisfying

where $H_{1}$ and $H_{2}$ are known constant matrices of appropriate dimensions.

$\tau (m)$ and $\delta (m)$ are time-varying delays satisfying

The nonlinear function $f_{i}(\cdot)$ satisfies Inequality (2.14) due to its monotonic increasing nature with saturation

Furthermore, this inequality can be reformulated in matrix form as follows:

where $K = {\rm diag}(k_{1}, k_{2}, \cdots k_{n}), \forall x\in R^{n}$.

Lemma 2.2 (^[24]). For any constant matrix $W = W^{\rm T} > 0$, scalar $r > 1$, there is

Lemma 2.3 (^[24]). For every pair of vectors $X, Y \in R^{n}$, where $H > 0$ denotes a positive definite matrix, the subsequent inequality is guaranteed to exist

3.   Passivity criteria

The general form of a linear matrix inequality (LMI) is as follows:

where $x = (x_{1}, x_{2}, \cdots, x_{n})$ is a vector of real decision variables, and $A_{0}, A_{1}, \cdots, A_{n}$ are given symmetric matrices (typically real-valued). Note that the symbol $\geq 0$ indicates that the matrix is positive semidefinite, which means all its eigenvalues are non-negative. The solution set of LMI is a convex set, meaning that its solution space has favorable mathematical properties. However, this also limits the range of feasible solutions. This convexity makes LMI very useful in optimization problems, but it increases its constraint. LMI is often used in system stability analysis, such as for verifying the stability of the system through Lyapunov functions. The application environment requires LMI to satisfy strict mathematical conditions to ensure the stability of the system.

In the subsequent discussions, an exploration will be conducted into a stability criterion pertaining to Eq (2.12) when subjected to Dirichlet boundary conditions.

Theorem 3.1 (^[24]). For given scalars $\delta$ and $\tau$ satisfying Inequality (2.13), and under the condition $\omega (k, i) = 0$, the trivial solution of Eq (2.12) without molecular Brownian motion under Dirichlet boundary conditions is passive if there exist scalars $\gamma > 0$, matrices $\Lambda _{h}^{\rm T} = \Lambda_{h} > 0\; (h = 1, 2, 3)$, $P_{h'}^{\rm T} = P_{h'} > 0$, $R_{h'}^{\rm T} = R_{h'} > 0$, $Q_{h'}^{\rm T} = Q_{h'} > 0\; (h' = 1, 2, \cdots, 4)$, $T_{1}^{\rm T} = T_{1} > 0$, $S_{1}^{\rm T} = S_{1} > 0$ and $S_{2}^{\rm T} = S_{2} > 0$, such that the following linear matrix inequalities hold:

where

$\Sigma_{1, 1} = 2A^{\rm T} P_{1} A-P_{1} +P_{3} +Q_{1} +2\delta_{\min}(A-I)^{\rm T} R_{1}(A-I)-\frac{R_{1}}{\delta_{\min}} \\ {\kern 45pt} +2(\delta_{\max} -\delta_{\min})(A-I)^{\rm T} R_{2}(A-I)+2\lambda_{1}S_{1} ,$

$\Sigma_{1, 3} = \frac{R_{1}}{\delta_{\min}}$,

$\Sigma_{1, 5} = A^{\rm T} P_{1} +\delta_{\min}(A-I)^{\rm T} R_{1} +(\delta_{\max}-\delta_{\min})(A-I)^{\rm T}R_{2} -\Lambda_{3}$,

$\Sigma_{1, 6} = A^{\rm T} P_{1}\theta +\delta_{\min}(A-I)^{\rm T} R_{1}\theta +(\delta_{\max}-\delta_{\min})(A-I)^{\rm T}R_{2}\theta$,

$\Sigma_{1, 7} = A^{\rm T} P_{1}\theta +\delta_{\min}(A-I)^{\rm T} R_{1}\theta +(\delta_{\max}-\delta_{\min})(A-I)^{\rm T}R_{2}\theta$,

$\Sigma_{2, 2} = 5D^{\rm T} P_{2} D-P_{3} +5\tau_{\min} D^{\rm T} R_{3} D+5(\tau_{\max} -\tau_{\min})D^{\rm T} R_{4} D$,

$\Sigma_{3, 3} = Q_{2} -Q_{1} -\frac{R_{1}}{\delta_{\min}} -\frac{R_{2} }{\delta_{\max} -\delta_{\min}}$,

$\Sigma_{3, 4} = \frac{R_{2}}{\delta_{\max}-\delta_{\min}}$,

$\Sigma_{4, 4} = -Q_{2} -\frac{R_{2} }{\delta_{\max} -\delta_{\min}}$,

$\Sigma_{5, 5} = 2P_{1} + 2\delta_{\min}R_{1} + 2(\delta_{\max} -\delta_{\min})R_{2}- \gamma\Lambda_{3}$,

$\Sigma_{5, 6} = P_{1}\theta +\delta_{\min}R_{1}\theta + (\delta_{\max} -\delta_{\min})R_{2}\theta$,

$\Sigma_{5, 7} = P_{1}\theta +\delta_{\min}R_{1}\theta + (\delta_{\max} -\delta_{\min})R_{2}\theta$,

$\Sigma_{6, 6} = 2\theta^{\rm T}P_{1}\theta + 2 \delta_{\min}\theta^{\rm T}R_{1}\theta + 2(\delta_{\max} -\delta_{\min})R_{2}\theta - S_{1}$,

$\Sigma_{6, 7} = \theta^{\rm T}P_{1}\theta + \delta_{\min}\theta^{\rm T}R_{1}\theta + (\delta_{\max} -\delta_{\min})R_{2}\theta$,

$\Sigma_{7, 7} = 2\theta^{\rm T}P_{1}\theta + 2 \delta_{\min}\theta^{\rm T}R_{1}\theta + 2(\delta_{\max} -\delta_{\min})R_{2}\theta - S_{1}$,

$\prod_{1, 1} = 2C^{\rm T} P_{2}C-P_{2} + P_{4} + Q_{3} +2\tau_{\min} (C-I)^{\rm T} R_{3} (C-I)- \frac{R_{3} }{\tau_{\min}} \\ {\kern 45pt} +2(\tau_{\max} -\tau _{\min})(C-I)^{\rm T} R_{4} (C-I)+2\lambda_{2} S_{2} ,$

$\prod_{1, 3} = \frac{R_{3}}{\delta_{\min}}$,

$\prod_{1, 5} = K\Lambda_{1}$,

$\prod_{1, 7} = C^{\rm T} P_{2} +\tau_{\min} (C-I)^{\rm T} R_{3} +(\tau_{\max} -\tau_{\min})(C-I)^{\rm T} R_{4} -\Lambda_{3}$,

$\prod_{1, 8} = C^{\rm T} P_{2}\mu +\tau_{\min} (C-I)^{\rm T} R_{3} \mu +(\tau_{\max} -\tau_{\min})(C-I)^{\rm T} R_{4}\mu$,

$\prod_{1, 9} = C^{\rm T} P_{2}\mu +\tau_{\min} (C-I)^{\rm T} R_{3} \mu +(\tau_{\max} -\tau_{\min})(C-I)^{\rm T} R_{4}\mu$,

$\prod_{2, 2} = -P_{4}$,

$\prod_{2, 6} = K\Lambda_{2}$ $\prod_{3, 3} = Q_{3} -Q_{4} -\frac{R_{3}}{\tau_{\min} } -\frac{R_{4}}{\tau_{\max} -\tau_{\min}}$,

$\prod_{3, 4} = \frac{R_{4}}{\delta_{\max}-\delta_{\min}}$,

$\prod_{4, 4} = -Q_{4} -\frac{R_{4}}{\tau_{\max} -\tau_{\min}}$,

$\prod_{5, 5} = T_{1} -2\Lambda_{1}$,

$\prod_{6, 6} = 5B^{\rm T} P_{1} B+5\delta_{\min} B^{\rm T} R_{1} B+5(\delta_{\max} -\delta_{\min})B^{\rm T} R_{2} B-T_{1} -2\Lambda_{2}$,

$\prod_{7, 7} = 2P_{2} +2\tau_{\min}R_{3} +2(\tau_{\max} -\tau_{\min})R_{4} - \gamma\Lambda_{3}$,

$\prod_{7, 8} = P_{2}\mu + \tau_{\min}R_{3}\mu + (\tau_{\max} -\tau_{\min})R_{4}\mu$,

$\prod_{7, 9} = P_{2}\mu + \tau_{\min}R_{3}\mu + (\tau_{\max} -\tau_{\min})R_{4}\mu$,

$\prod_{8, 8} = 2\mu ^{\rm T}P_{2}\mu + 2\tau_{\min}\mu ^{\rm T}R_{3}\mu + 2(\delta_{\max} -\delta_{\min})\mu ^{\rm T}R_{4}\mu - S_{2}$,

$\prod_{8, 9} = \mu ^{\rm T}P_{2}\mu + \tau_{\min}\mu ^{\rm T}R_{3}\mu + (\delta_{\max} -\delta_{\min})\mu ^{\rm T}R_{4}\mu$,

$\prod_{9, 9} = 2\mu ^{\rm T}P_{2}\mu + 2\tau_{\min}\mu ^{\rm T}R_{3}\mu + 2(\delta_{\max} -\delta_{\min})\mu ^{\rm T}R_{4}\mu - S_{2}$.

Proof. Set

Consider a Lyapunov-Krasovskii functional candidate defined as

where

Define $\Delta V(k, i) = V(k+1, i)-V(k, i)$, then

where

Taking Eq (3.4) into consideration, for diagonal matrices $\Lambda_{1} > 0$, $\Lambda_{2} > 0$, there is

Using Lemma 2.2, we have

The following expression is given to show the passivity property

and it is obvious that

where $k\in N$, such that for the original time $k_{0} = 1$, set the initial states of $m$ and $p$ as $m(1, i) = p(1, i) = 0$, so that taking (3.6)–(3.10) into account, we can obtain

Considering the spatial continuity and limitations of GRNs, we can obtain

where

Then, by combining Eq (3.11)–(3.27), we can obtain

where $\xi = \left[\begin{aligned} & m^{\rm T}(k, i), m^{\rm T}(k-\delta (k), i), m^{\rm T} (k-\delta_{\min}, i), m^{\rm T}(k-\delta_{\max}, i), p^{\rm T}(k, i), \\ & p^{\rm T}(k-\tau(k), i), p^{\rm T}(k-\tau_{\min}, i), p^{\rm T}(k-\tau_{\max}, i), f^{\rm T}(p(k, i)), \\ & f^{\rm T}(p(k-\tau(k), i)), u^{\rm T}(k, i), v^{\rm T}(k, i), m^{\rm T}(k, i+1), m^{\rm T}(k, i-1), \\ & p^{\rm T} (k, i+1), p^{\rm T}(k, i-1) \end{aligned} \right].$

This completes the proof. □

Theorem 3.2. For given scalars $\delta$ and $\tau$ satisfying inequality (2.13), the trivial solution of Eq (2.12) with molecular Brownian motion under Dirichlet boundary conditions is passive if there exist scalars $\gamma > 0$ and $\rho > 0$, matrices $\Lambda_{h}^{\rm T} = \Lambda_{h} > 0\; (h = 1, 2)$, $P_{h'}^{\rm T} = P_{h'} > 0$, $R_{h'}^{\rm T} = R_{h'} > 0$ and $Q_{h'}^{\rm T} = Q_{h'} > 0\; (h' = 1, \cdots, 4)$, $T_{1}^{\rm T} = T_{1} > 0$, $S_{1}^{\rm T} = S_{1} > 0$ and $S_{2}^{\rm T} = S_{2} > 0$, such that the following LMIs hold:

where

$\sum_{11} = A^{\rm T} P_{1} A+\delta_{\min } (A-I)^{\rm T} R_{1} (A-I)+(\delta_{\max} -\delta_{\min})(A-I)^{\rm T} R_{2} (A-I) \\ {\kern 45pt} -P_{1} +P_{3} +Q_{1} -R_{1} /\delta_{\min} +2\lambda_{1} S_{1},$

$\sum_{13} = \frac{R_{1}}{\delta_{\min}}$,

$\sum_{15} = A^{\rm T} P_{1} B+\delta_{\min} (A-I)^{\rm T} R_{1} B+(\delta_{\max} -\delta_{\min})(A-I)^{\rm T} R_{2} B$,

$\sum_{16} = A^{\rm T} P_{1} +\delta_{\min} (A-I)^{\rm T} R_{1} +(\delta_{\max} -\delta_{\min})(A-I)^{\rm T} R_{2} +\Lambda_{3}$,

$\sum_{17} = A^{\rm T} P_{1} \theta +\delta_{\min} (A-I)^{\rm T} R_{1} \theta +(\delta_{\max} -\delta_{\min})(A-I)^{\rm T} R_{2} \theta$,

$\sum_{18} = A^{\rm T} P_{1} \theta +\delta_{\min} (A-I)^{\rm T} R_{1} \theta +(\delta_{\max} -\delta_{\min})(A-I)^{\rm T} R_{2} \theta$,

$\sum_{22} = 3D^{\rm T} P_{2} D-P_{3} +3\tau_{\min} D^{\rm T} R_{3} D+3(\tau_{\max} -\tau_{\min})D^{\rm T} R_{4} D$,

$\sum_{33} = Q_{2} -Q_{1} -\frac{R_{1}}{\delta_{\min}} -\frac{R_{2}}{\tau_{\max} -\tau_{\min}}$,

$\sum_{34} = \frac{R_{2}}{\delta_{\max} -\delta_{\min}}$,

$\sum_{44} = -Q_{2} -\frac{R_{2}}{\delta_{\max} -\delta_{\min}}$,

$\sum_{55} = B^{\rm T} P_{1} B+\delta_{\min} B^{\rm T} R_{1} B+(\delta_{\max} -\delta_{\min})B^{\rm T} R_{2} B-T_{1} +\Lambda_{2}$,

$\sum_{56} = B^{\rm T} P_{1} +\delta_{\min} B^{\rm T} R_{1} +(\delta _{\max } -\delta_{\min})B^{\rm T} R_{2}$,

$\sum_{57} = B^{\rm T} P\theta +\delta_{\min} B^{\rm T} R_{1} \theta +(\delta_{\max} -\delta_{\min})B^{\rm T} R_{2} \theta$,

$\sum_{58} = B^{\rm T} P\theta +\delta_{\min} B^{\rm T} R_{1} \theta +(\delta_{\max} -\delta_{\min})B^{\rm T} R_{2} \theta$,

$\sum_{66} = P_{1} +\delta_{\min} R_{1} +(\delta_{\max} -\delta_{\min})R_{2} -\gamma \Lambda_{3}$,

$\sum_{67} = P_{1} \theta +\delta_{\min } R_{1} \theta +(\delta_{\max } -\delta _{\min})R_{2} \theta$,

$\sum_{68} = P_{1} \theta +\delta_{\min } R_{1} \theta +(\delta_{\max } -\delta _{\min})R_{2} \theta$,

$\sum_{77} = \theta^{\rm T} P_{1} \theta +\delta_{\min} \theta^{\rm T} R_{1} \theta +(\delta_{\max} -\delta_{\min})\theta^{\rm T} R_{2} \theta -S_{1}$,

$\sum_{78} = \theta^{\rm T} P_{1} \theta +\delta_{\min} \theta^{\rm T} R_{1} \theta +(\delta_{\max} -\delta_{\min})\theta^{\rm T} R_{2} \theta$,

$\sum_{88} = \theta^{\rm T} P_{1} \theta +\delta_{\min} \theta^{\rm T} R_{1} \theta +(\delta_{\max} -\delta_{\min})\theta^{\rm T} R_{2} \theta -S_{1}$,

$\Pi_{11} = 2C^{\rm T} P_{2} C-P_{2} + P_{4} +Q_{3} +2\tau_{\min} (C-I)^{\rm T} R_{3} (C-I)-R_{3} /\tau_{\min} \\ {\kern 45pt}+2(\tau_{\max} -\tau_{\min})(C-I)^{\rm T} R_{4} (C-I)+\rho H_{1} +2\lambda_{2} S_{2} ,$

$\prod_{13} = R_{3} /\tau_{\min}$,

$\prod_{15} = K \Lambda_{1}$,

$\prod_{16} = C^{\rm T} P_{2} +\tau_{\min} (C-I)^{\rm T} R_{3} +(\tau_{\max} -\tau_{\min})(C-I)^{\rm T} R_{4} +\Lambda_{3}$,

$\prod_{17} = C^{\rm T} P_{2} \mu_{1} +\tau_{\min} (C-I)^{\rm T} R_{3} \mu_{1} +(\tau_{\max} -\tau_{\min})(C-I)^{\rm T} R_{4} \mu_{1}$,

$\prod_{18} = C^{\rm T} P_{2} \mu_{2} +\tau_{\min} (C-I)^{\rm T} R_{3} \mu_{2} +(\tau_{\max} -\tau_{\min})(C-I)^{\rm T} R_{4} \mu_{2}$,

$\prod_{22} = -P_{4} +K^{\rm T} \Lambda_{2} K+\rho H_{2}$,

$\prod_{33} = Q_{3} -Q_{4} -R_{3} /\tau_{\min} -R_{4} /(\tau_{\max} -\tau_{\min})$,

$\prod_{34} = R_{4} /(\tau_{\max} -\tau_{\max})$,

$\prod_{44} = -Q_{4} -R_{4} /(\tau_{\max} -\tau_{\min})$,

$\prod_{55} = -2\Lambda_{1} +T_{1}$,

$\prod_{66} = P_{2} +\tau_{\max} R_{3} +(\tau_{\max} -\tau_{\min})R_{4} -\gamma \Lambda_{3}$,

$\prod_{67} = P_{2} \mu +\tau_{\min } R_{3} \mu +(\tau_{\max } -\tau_{\min })R_{4} \mu$,

$\prod_{68} = P_{2} \mu +\tau_{\min } R_{3} \mu +(\tau_{\max } -\tau_{\min })R_{4} \mu$,

$\prod_{77} = \mu^{\rm T} P_{2} \mu +\tau_{\min} \mu^{\rm T} R_{3} \mu +(\tau_{\max} -\tau_{\min})\mu^{\rm T} R_{4} \mu -S_{2}$,

$\prod_{78} = \mu^{\rm T} P_{2} \mu +\tau_{\min} \mu^{\rm T} R_{3} \mu +(\tau_{\max} -\tau_{\min})\mu^{\rm T} R_{4} \mu$,

$\prod_{88} = \mu^{\rm T} P_{2} \mu +\tau_{\min} \mu^{\rm T} R_{3} \mu +(\tau_{\max} -\tau_{\min})\mu^{\rm T} R_{4} \mu-S_{2}$.

Proof. Using Lemma 2.3, we have

Noting Inequality (2.13), it can be seen that

Then combining Eqs (3.11) to (3.25) with Eqs (3.32) and (3.33), we get

where

$\xi_{1} = \left[\begin{aligned} & m^{\rm T} (k, i), m^{\rm T}\left(k-\delta (k), i\right), m^{\rm T}(k-\delta_{\min}, i), m^{\rm T}(k-\delta_{\max}, i), \\ & f^{\rm T}\left(p\left(k-\tau (k), i\right)\right), u^{\rm T}(k, i), m^{\rm T}(k, i+1), m^{\rm T}(k, i-1) \end{aligned} \right],$

$\xi_{2} = \left[\begin{aligned} & p^{\rm T}(k, i), p^{\rm T}\left(k-\tau (k), i\right), p^{\rm T}(k-\tau_{\min}, i), p^{\rm T}(k-\tau_{\max}, i), \\ & f^{\rm T}\left(p(k, i)\right), v^{\rm T}(k, i), p^{\rm T}(k, i+1), p^{\rm T}(k, i-1) \end{aligned} \right].$

The proof of the theorem is now complete. □

Remark 3.3. In real genetic regulatory networks, the functions and operations of mRNA and proteins are often affected by noise, which includes physical noise and biological noise. Among these, Brownian motion represents a form of physical noise that arises from the random movements of large molecules such as mRNA and proteins within the cell. This randomness contributes to fluctuations in molecular concentrations and interactions, thereby affecting the overall functionality and stability of the network. Theorem 3.1 presents a foundational understanding of these networks without accounting for certain types of noise. In contrast, Theorem 3.2 takes into account the impact of physical noise generated by Brownian motion on the passivity of the model to a greater extent. This not only enhances the model's predictive power but also provides deeper insights into the role of physical noise in shaping network behavior.

4.   Simulation results

In this section, we demonstrate the validity of Theorems 3.1 and 3.2 through two examples. Example 4.1 showcases the passivity of the proposed GRNs model without Brownian motion, while Example 4.2 illustrates the passivity of the proposed GRNs model with Brownian motion included.

Example 4.1. In this section, we consider a GRN (2.12) with $5$ nodes, $L = 100$. The parameters are assumed to be

$A = C = {\rm diag}\left(0.1, 0.1, 0.1, 0.1, 0.1\right), \; D = \left(0.8, 0.8, 0.8, 0.8, 0.8\right)$,

$B = \left[\begin{array}{ccccc} {0} & {-0.5} & {0.5} & {0} & {0} \\ {-0.5} & {0} & {0} & {0.5} & {0.5} \\ {0} & {0.5} & {0} & {0} & {0} \\ {0.5} & {-0.5} & {0} & {0} & {0} \\ {0} & {0} & {0} & {0.5} & {0} \end{array} \right], \; K = \left[\begin{array}{ccccc} {0.65} & {0} & {0} & {0} & {0} \\ {0} & {0.65} & {0} & {0} & {0} \\ {0} & {0} & {0.65} & {0} & {0} \\ {0} & {0} & {0} & {0.65} & {0} \\ {0} & {0} & {0} & {0} & {0.65} \end{array} \right]$.

$u_{1}(k, i) = 0.5\left(1-\frac{i-\frac{L}{2} }{\frac{L}{2}}\right)^{2} \sin(2k)\times m_{1}(k, i)$, $u_{2} (k, i) = 0.5\left(1.2-\frac{i-\frac{L}{2}}{\frac{L}{2}} \right)^{2} \cos(2k)\times m_{2}(k, i)$,

$u_{3} (k, i) = 0.5\left(0.5-\frac{i-\frac{L}{2}}{\frac{L}{2}} \right)^{2} \sin(k)\times m_{3}(k, i)$, $u_{4} (k, i) = 0.5\left(0.2-\frac{i-\frac{L}{2} }{\frac{L}{2}} \right)^{2} \cos(0.2k)\times m_{4}(k, i)$,

$u_{5} (k, i) = 0.5\left(0.5-\frac{i-\frac{L}{2}}{\frac{L}{2}} \right)^{2} \sin(0.3k)\times m_{5}(k, i)$, $v_{1}(k, i) = 0.5\left(1.2-\frac{i-\frac{L}{2}}{\frac{L}{2}} \right)^{2} \sin(3k)\times p_{1}(k, i)$,

$v_{2}(k, i) = 0.5\left(1-\frac{i-\frac{L}{2}}{\frac{L}{2}} \right)^{2} \cos(k)\times p_{2}(k, i)$, $v_{3}(k, i) = 0.5\left(0.7-\frac{i-\frac{L}{2}}{\frac{L}{2}} \right)^{2} \cos(2k)\times p_{3}(k, i)$,

$v_{4}(k, i) = 0.5\left(1-\frac{i-\frac{L}{2}}{\frac{L}{2}} \right)^{2} \cos(k)\times p_{4}(k, i)$, $v_{5}(k, i) = 0.5\left(0.7-\frac{i-\frac{L}{2}}{\frac{L}{2}} \right)^{2} \cos(2k)\times p_{5}(k, i)$,

where $f(x) = x^{2} /(1+x^{2})$, the time delays $\delta (k) = 4+2\sin (k\pi /2)$ and $\tau(k) = 4+\sin (k\pi /2)$, so that $\delta _{\min } = 2$, $\delta_{\max } = 6$, $\tau_{\min} = 3$, $\tau_{\max} = 5$, coupling coefficient $\theta = \mu = 0.2$.

The simulation result of trajectories of mRNA and protein concentrations for Example 4.1 are shown in Figures 2–6.

If $\phi'_{1} = \phi'_{2} = \phi'\in \left\{0.1, 0.3, 0.5, 0.7, 0.9, 1.0 \right\}$, $\phi_{1} = \phi_{2} = \phi > 1$. The maximum delay $\tau_{{\rm max}} = \delta_{\max}$ is obtained based on (^[21], Theorem 1), (^[25], Remark 2), (^[22], Theorem 1), and Theorem 3.1; the following illustrations are made for Table 1:

● If $\mu \le 0.7$, the LMI conditions outlined in (^[21], Theorem 1), (^[25], Remark 2), (^[22], Theorem 1), and Theorem 3.1 are considered to be achievable.

● When $\mu = 0.9$, the LMI conditions given in (^[22], Theorem 1) are feasible, but the LMI conditions presented in (^[21], Theorem 1), (^[25], Remark 2), and Theorem 3.1 are $\mathrm{\infty}$.

● When $\mu = 1$, the LMI conditions discussed in (^[21], Theorem 1) and (^[25], Remark 2) are feasible, the LMI conditions discussed in (^[22], Theorem 1) are unfeasible, and only the LMI conditions discussed in Theorem 3.1 are $\mathrm{\infty}$.

Therefore, within the range of $\mu \le 1$, Theorem 3.1 of this paper exhibits a reduced level of conservativeness compared to (^[21], Theorem 1), (^[25], Remark 2), and (^[22], Theorem 1).

Example 4.2. Considering the Brownian motion for GRNs (2.12), $H_{1} = H_{2} = 1$, $\lambda _{1} = \lambda _{2} = 2$, by using the Toolbox YALMIP in MATLAB to solve the LMIs (3.2) and (3.3), we can obtain the following feasible solution:

Therefore, it can be concluded by Examples 4.1 and 4.2 that Eq (2.8) is stochastically passive.

Figure 7 presents the boundary between the passive and non-passive behavior regions of the proposed GRN model at $\mu = 0.55$ and $\theta = 0.60$. It can be clearly observed that when the values of the coupling coefficients $\theta$ and $\mu$ fall within the parameter range defined by Region Ⅰ, specifically $0.05 \leq \mu \leq 0.55$ and $0.05 \leq \theta \leq 0.60$, the proposed discrete-time gene regulatory network model exhibits passive behavior. Conversely, when the values of $\theta$ and $\mu$ fall within the range defined by Region Ⅱ, the proposed GRN model does not exhibit passive behavior. The dashed line in Figure 7 precisely delineates the boundary between the regions of passivity and non-passivity, distinguishing Region Ⅰ from Region Ⅱ. This boundary was determined through extensive computer software simulations, identifying the critical values of $\theta$ and $\mu$ at which the system transitions from passive to non-passive behavior. The shapes of the regional boundaries reflect the complex interplay between the model parameters and their impact on the system's passivity properties.

Remark 4.3. In real biological environments, the concentrations of mRNA and proteins are indeed influenced by a wide array of factors, including gene mutations, interactions between cells, biological enzymes, hormone levels, and environmental conditions. These factors contribute to the complex dynamics and heterogeneity observed in biological systems. However, our proposed model focuses on a more simplified mathematical representation of these processes, specifically considering the interactions between mRNA and proteins, the effects of spatial diffusion terms, and the impact of time delays on their concentrations. While this mathematical simplification allows us to derive analytical results and gain insights into the passivity and stability of the system, it also introduces limitations in the model's ability to fully capture the complexity of real biological environments. For instance, our model does not account for the potential effects of gene mutations or the intricate interactions between cells in a tissue or organ. Additionally, we do not consider the role of biological enzymes or hormone levels, which can significantly impact gene expression and protein synthesis.

It is important to recognize these limitations and to understand that our model represents a first step toward understanding the dynamics of gene regulatory networks in a mathematical framework. Future work could build upon our findings by incorporating additional factors that contribute to the complexity of real biological systems. Such extensions could provide even deeper insights into the mechanisms underlying gene expression and protein synthesis, as well as their implications for health and disease.

5.   Conclusions

In this groundbreaking study, we propose a pioneering approach to reaction-diffusion mechanisms into discrete-time GRNs. Our primary focus is to explore the robust passivity of these networks under the influence of time-varying delay and Dirichlet boundary conditions using advanced Lyapunov-Krasovskii functions. Moreover, we delve into the investigation of asymptotic passivity for GRNs that incorporate reaction-diffusion terms along with Brownian motion. This unique exploration is unprecedented and marks the first-ever attempt to introduce reaction-diffusion into discrete-time GRNs. To demonstrate the efficacy and validity of our novel method, we present comprehensive numerical examples and simulation results. These findings showcase the correctness and effectiveness of our approach, solidifying its potential impact on the field of gene regulatory network modeling and analysis. In the future, the asynchronous behavior of discrete-time systems with reaction-diffusion coupling could be studied, in order to understand how different update schedules affect the passivity and stability of the network. In addition, we can also try to extend the proposed model to more complex spatial domains, such as those with irregular geometries or anisotropic diffusion, which will provide further insights into the spatial organization of gene regulatory networks.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was funded by Post-doctoral Natural Science Foundation of China under grant number 2022M720650.

Conflict of interest

The authors declare there is no conflicts of interest.