Models Of Regulatory Networks Research Articles

On the evolutionary origin of discrete phenotypic plasticity.

Phenotypic plasticity provides an attractive strategy for adapting to various environments, but the evolutionary mechanism of the underlying genetic system is poorly understood. We use a simple gene regulatory network model to explore how a species acquires phenotypic plasticity, particularly focusing on discrete phenotypic plasticity, which has been difficult to explain by quantitative genetic models. Our approach employs a population genetic framework that integrates the developmental process, where each individual undergoes growth to develop its phenotype, which subsequently becomes subject to selection pressures. Our model considers two alternative types of environments, with the gene regulatory network including a sensor gene that turns on and off depending on the type of environment. With this assumption, we demonstrate that the system gradually adapts by acquiring the ability to produce two distinct optimum phenotypes under two types of environments without changing genotype, resulting in phenotypic plasticity. We find that the resulting plasticity is often discrete after a lengthy period of evolution. Our results suggest that gene regulatory networks have a notable capacity to flexibly produce various phenotypes in response to environmental changes. This study also shows that the evolutionary dynamics of phenotype may differ significantly between mechanistic-based developmental models and quantitative genetics models, suggesting the utility of incorporating gene regulatory networks into evolutionary models.

Uniform convergence of solutions to stochastic hybrid models of gene regulatory networks

In a recent paper by Kurasov et al. (Math Biosci 305:170–177, 2018), a hybrid gene regulatory network was proposed to model gene expression dynamics by using a stochastic system of coupled partial differential equations. In more recent work, the existence and strong convergence of the solutions to equilibrium were proven. In this article, we improve upon their result by showing that the convergence rate is independent of the initial state, therefore proving that the solutions converge not only strongly but even uniformly to equilibrium. To this end, we make use of a recent convergence theorem for stochastic, irreducible semigroups that contain partial integral operators.

D-SPIN constructs gene regulatory network models from multiplexed scRNA-seq data revealing organizing principles of cellular perturbation response.

Gene regulatory networks within cells modulate the expression of the genome in response to signals and changing environmental conditions. Reconstructions of gene regulatory networks can reveal the information processing and control principles used by cells to maintain homeostasis and execute cell-state transitions. Here, we introduce a computational framework, D-SPIN, that generates quantitative models of gene regulatory networks from single-cell mRNA-seq datasets collected across thousands of distinct perturbation conditions. D-SPIN models the cell as a collection of interacting gene-expression programs, and constructs a probabilistic model to infer regulatory interactions between gene-expression programs and external perturbations. Using large Perturb-seq and drug-response datasets, we demonstrate that D-SPIN models reveal the organization of cellular pathways, sub-functions of macromolecular complexes, and the logic of cellular regulation of transcription, translation, metabolism, and protein degradation in response to gene knockdown perturbations. D-SPIN can also be applied to dissect drug response mechanisms in heterogeneous cell populations, elucidating how combinations of immunomodulatory drugs can induce novel cell states through additive recruitment of gene expression programs. D-SPIN provides a computational framework for constructing interpretable models of gene-regulatory networks to reveal principles of cellular information processing and physiological control.

Identification of Potential lncRNA-miRNA-mRNA Regulatory Network Contributing to Arrhythmogenic Right Ventricular Cardiomyopathy.

Arrhythmogenic right ventricular cardiomyopathy (ARVC) can lead to sudden cardiac death and life-threatening heart failure. Due to its high fatality rate and limited therapies, the pathogenesis and diagnosis biomarker of ARVC needs to be explored urgently. This study aimed to explore the lncRNA-miRNA-mRNA competitive endogenous RNA (ceRNA) network in ARVC. The mRNA and lncRNA expression datasets obtained from the Gene Expression Omnibus (GEO) database were used to analyze differentially expressed mRNA (DEM) and lncRNA (DElnc) between ARVC and non-failing controls. Differentially expressed miRNAs (DEmiRs) were obtained from the previous profiling work. Using starBase to predict targets of DEmiRs and intersecting with DEM and DElnc, a ceRNA network of lncRNA-miRNA-mRNA was constructed. The DEM and DElnc were validated by real-time quantitative PCR in human heart tissue. Protein-protein interaction network and weighted gene co-expression network analyses were used to identify hub genes. A logistic regression model for ARVC diagnostic prediction was established with the hub genes and their ceRNA pairs in the network. A total of 448 DEMs (282 upregulated and 166 downregulated) were identified, mainly enriched in extracellular matrix and fibrosis-related GO terms and KEGG pathways, such as extracellular matrix organization and collagen fibril organization. Four mRNAs and two lncRNAs, including COL1A1, COL5A1, FBN1, BGN, XIST, and LINC00173 identified through the ceRNA network, were validated by real-time quantitative PCR in human heart tissue and used to construct a logistic regression model. Good ARVC diagnostic prediction performance for the model was shown in both the training set and the validation set. The potential lncRNA-miRNA-mRNA regulatory network and logistic regression model established in our study may provide promising diagnostic methods for ARVC.

Flexible modeling of regulatory networks improves transcription factor activity estimation

Transcriptional regulation plays a crucial role in determining cell fate and disease, yet inferring the key regulators from gene expression data remains a significant challenge. Existing methods for estimating transcription factor (TF) activity often rely on static TF-gene interaction databases and cannot adapt to changes in regulatory mechanisms across different cell types and disease conditions. Here, we present a new algorithm - Transcriptional Inference using Gene Expression and Regulatory data (TIGER) - that overcomes these limitations by flexibly modeling activation and inhibition events, up-weighting essential edges, shrinking irrelevant edges towards zero through a sparse Bayesian prior, and simultaneously estimating both TF activity levels and changes in the underlying regulatory network. When applied to yeast and cancer TF knock-out datasets, TIGER outperforms comparable methods in terms of prediction accuracy. Moreover, our application of TIGER to tissue- and cell-type-specific RNA-seq data demonstrates its ability to uncover differences in regulatory mechanisms. Collectively, our findings highlight the utility of modeling context-specific regulation when inferring transcription factor activities.

Genome-wide analysis of long noncoding RNAs as cis-acting regulators of transcription factor-encoding genes in IgA nephropathy.

IgA nephropathy (IgAN) is the most common form of primary glomerulonephritis in the world, but the disease pathogenesis noncoding is yet to be elucidated. Previous studies have revealed regulatory functions for long noncoding RNA (lncRNA) in various diseases; however, the roles of lncRNA in IgAN and regulation of transcription factors (TFs) have been scarcely investigated. Renal tissue samples (n = 5) from patients with IgAN and control samples (n = 4) were collected and RNA sequencing (RNA-seq) was performed. Four software programs were employed for lncRNA prediction. GO (Gene Ontology)/KEGG (Kyoto Encyclopedia of Genes and Genomes) were employed for analysis of the identified differentially expressed genes (DEGs). A regulatory network model of DE lncRNA-TF-DEG was developed, and the levels of expression of key lncRNAs, TFs, and corresponding target genes were assessed using qRT-PCR and immunofluorescence. The current study identified 674 upregulated and 1,011 downregulated DE mRNAs and 260 upregulated and 232 downregulated DE lncRNAs in IgAN samples compared with control samples. The upregulated DE mRNAs showed enrichment in cell adhesion and collagen glial fiber organization pathways. The DE lncRNAs-DE mRNAs showing co-expression are associated with transmembrane transport. A novel regulatory network model of lncRNA-TF-DEG has been developed. This study identified seven TFs that are cis-regulated by 6 DE lncRNAs, and show co-expression with 132 DEGs (correlation coefficient ≥ 0.8, P ≤ 0.01), generating 158 pairs that showed co-expression. The lncRNAs NQO1-DT and RP5-1057120.6 were found to be highly expressed in IgAN samples. The TFs vitamin D Receptor (VDR) and NFAT5, along with their target genes were also aberrantly expressed. Key lncRNAs and TFs centrally associated with IgAN have been identified in this study. A regulatory network model of lncRNA-TF-mRNA was constructed. Further studies on the genes identified herewith could provide insight into the pathogenesis of IgAN.

Characterization of clinical data for patient stratification in moderate osteoarthritis with support vector machines, regulatory network models, and verification against osteoarthritis Initiative data

Knee osteoarthritis (OA) diagnosis is based on symptoms, assessed through questionnaires such as the WOMAC. However, the inconsistency of pain recording and the discrepancy between joint phenotype and symptoms highlight the need for objective biomarkers in knee OA diagnosis. To this end, we study relationships among clinical and molecular data in a cohort of women (n = 51) with Kellgren–Lawrence grade 2–3 knee OA through a Support Vector Machine (SVM) and a regulation network model. Clinical descriptors (i.e., pain catastrophism, depression, functionality, joint pain, rigidity, sensitization and synovitis) are used to classify patients. A Youden’s test is performed for each classifier to determine optimal binarization thresholds for the descriptors. Thresholds are tested against patient stratification according to baseline WOMAC data from the Osteoarthritis Initiative, and the mean accuracy is 0.97. For our cohort, the data used as SVM inputs are knee OA descriptors, synovial fluid proteomic measurements (n = 25), and transcription factor activation obtained from regulatory network model stimulated with the synovial fluid measurements. The relative weights after classification reflect input importance. The performance of each classifier is evaluated through ROC-AUC analysis. The best classifier with clinical data is pain catastrophism (AUC = 0.9), highly influenced by funcionality and pain sensetization, suggesting that kinesophobia is involved in pain perception. With synovial fluid proteins used as input, leptin strongly influences every classifier, suggesting the importance of low-grade inflammation. When transcription factors are used, the mean AUC is limited to 0.608, which can be related to the pleomorphic behaviour of osteoarthritic chondrocytes. Nevertheless, funcionality has an AUC of 0.7 with a decisive importance of FOXO downregulation. Though larger and longitudinal cohorts are needed, this unique combination of SVM and regulatory network model shall help to stratify knee OA patients more objectively.

Biologically informed NeuralODEs for genome-wide regulatory dynamics

BackgroundGene regulatory network (GRN) models that are formulated as ordinary differential equations (ODEs) can accurately explain temporal gene expression patterns and promise to yield new insights into important cellular processes, disease progression, and intervention design. Learning such gene regulatory ODEs is challenging, since we want to predict the evolution of gene expression in a way that accurately encodes the underlying GRN governing the dynamics and the nonlinear functional relationships between genes. Most widely used ODE estimation methods either impose too many parametric restrictions or are not guided by meaningful biological insights, both of which impede either scalability, explainability, or both.ResultsWe developed PHOENIX, a modeling framework based on neural ordinary differential equations (NeuralODEs) and Hill-Langmuir kinetics, that overcomes limitations of other methods by flexibly incorporating prior domain knowledge and biological constraints to promote sparse, biologically interpretable representations of GRN ODEs. We tested the accuracy of PHOENIX in a series of in silico experiments, benchmarking it against several currently used tools. We demonstrated PHOENIX’s flexibility by modeling regulation of oscillating expression profiles obtained from synchronized yeast cells. We also assessed the scalability of PHOENIX by modeling genome-scale GRNs for breast cancer samples ordered in pseudotime and for B cells treated with Rituximab.ConclusionsPHOENIX uses a combination of user-defined prior knowledge and functional forms from systems biology to encode biological “first principles” as soft constraints on the GRN allowing us to predict subsequent gene expression patterns in a biologically explainable manner.

Temporal multiscale modeling of biochemical regulatory networks: Calcium-regulated hepatocyte lipid and glucose metabolism

Hepatocyte lipid and glucose metabolism is regulated not only by major hormones like insulin and glucagon but also by many other factors, including calcium ions. Recently, mitochondria-associated membrane (MAM) dysfunction combined with incorrect IP3-receptor regulation has been shown to result in abnormal calcium signaling in hepatocytes. This dysfunction could further lead to hepatic metabolism pathology. However, the exact contribution of MAM dysfunction, incorrect IP3-receptor regulation and insulin resistance to the calcium-insulin-glucagon interplay is not understood yet. In this work, we analyze the role of abnormal calcium signaling and insulin dysfunction in hepatocytes by proposing a model of hepatocyte metabolic regulatory network with a detailed focus on the model construction details besides the biological aspect. In this work, we analyze the role of abnormal calcium signaling and insulin dysfunction in hepatocytes by proposing a model of hepatocyte metabolic regulatory network. We focus on the model construction details, model validation, and predictions. We describe the dynamic regulation of signaling processes by sigmoid Hill function. In particular, we study the effect of both the Hill function slope and the distance between Hill function extremes on metabolic processes in hepatocytes as a model of nonspecific insulin dysfunction. We also address the significant time difference between characteristic time of glucose hepatic processing and a typical calcium oscillation period in hepatocytes. Our modeling results show that calcium signaling dysfunction results in an abnormal increase in postprandial glucose levels, an abnormal glucose decrease in fasting, and a decreased amount of stored glycogen. An insulin dysfunction of glucose phosphorylation, glucose dephosphorylation, and glycogen breakdown also cause a noticeable effect. We also get some insight into the so-called hepatic insulin resistance paradox, confirming the hypothesis regarding indirect insulin action on hepatocytes via dysfunctional adipocyte lipolysis.

Effects of Small Random Perturbations in the Extended Glass–Kauffman Model of Gene Regulatory Networks

Effects of Small Random Perturbations in the Extended Glass–Kauffman Model of Gene Regulatory Networks

Data-driven modeling of core gene regulatory network underlying leukemogenesis in IDH mutant AML

Acute myeloid leukemia (AML) is characterized by uncontrolled proliferation of poorly differentiated myeloid cells, with a heterogenous mutational landscape. Mutations in IDH1 and IDH2 are found in 20% of the AML cases. Although much effort has been made to identify genes associated with leukemogenesis, the regulatory mechanism of AML state transition is still not fully understood. To alleviate this issue, here we develop a new computational approach that integrates genomic data from diverse sources, including gene expression and ATAC-seq datasets, curated gene regulatory interaction databases, and mathematical modeling to establish models of context-specific core gene regulatory networks (GRNs) for a mechanistic understanding of tumorigenesis of AML with IDH mutations. The approach adopts a new optimization procedure to identify the top network according to its accuracy in capturing gene expression states and its flexibility to allow sufficient control of state transitions. From GRN modeling, we identify key regulators associated with the function of IDH mutations, such as DNA methyltransferase DNMT1, and network destabilizers, such as E2F1. The constructed core regulatory network and outcomes of in-silico network perturbations are supported by survival data from AML patients. We expect that the combined bioinformatics and systems-biology modeling approach will be generally applicable to elucidate the gene regulation of disease progression.

Effects of mRNA conformational switching on translational noise in gene circuits.

Intragenic translational heterogeneity describes the variation in translation at the level of transcripts for an individual gene. A factor that contributes to this source of variation is the mRNA structure. Both the composition of the thermodynamic ensemble, i.e., the stationary distribution of mRNA structures, and the switching dynamics between those play a role. The effect of the switching dynamics on intragenic translational heterogeneity remains poorly understood. We present a stochastic translation model that accounts for mRNA structure switching and is derived from a Markov model via approximate stochastic filtering. We assess the approximation on various timescales and provide a method to quantify how mRNA structure dynamics contributes to translational heterogeneity. With our approach, we allow quantitative information on mRNA switching from biophysical experiments or coarse-grain molecular dynamics simulations of mRNA structures to be included in gene regulatory chemical reaction network models without an increase in the number of species. Thereby, our model bridges a gap between mRNA structure kinetics and gene expression models, which we hope will further improve our understanding of gene regulatory networks and facilitate genetic circuit design.

Biologically meaningful regulatory logic enhances the convergence rate in Boolean networks and bushiness of their state transition graph.

Boolean models of gene regulatory networks (GRNs) have gained widespread traction as they can easily recapitulate cellular phenotypes via their attractor states. Their overall dynamics are embodied in a state transition graph (STG). Indeed, two Boolean networks (BNs) with the same network structure and attractors can have drastically different STGs depending on the type of Boolean functions (BFs) employed. Our objective here is to systematically delineate the effects of different classes of BFs on the structural features of the STG of reconstructed Boolean GRNs while keeping network structure and biological attractors fixed, and explore the characteristics of BFs that drive those features. Using $10$ reconstructed Boolean GRNs, we generate ensembles that differ in BFs and compute from their STGs the dynamics' rate of contraction or 'bushiness' and rate of 'convergence', quantified with measures inspired from cellular automata (CA) that are based on the garden-of-Eden (GoE) states. We find that biologically meaningful BFs lead to higher STG 'bushiness' and 'convergence' than random ones. Obtaining such 'global' measures gets computationally expensive with larger network sizes, stressing the need for feasible proxies. So we adapt Wuensche's $Z$-parameter in CA to BFs in BNs and provide four natural variants, which, along with the average sensitivity of BFs computed at the network level, comprise our descriptors of local dynamics and we find some of them to be good proxies for bushiness. Finally, we provide an excellent proxy for the 'convergence' based on computing transient lengths originating at random states rather than GoE states.

Plasticity-led and mutation-led evolutions are different modes of the same developmental gene regulatory network.

The standard theory of evolution proposes that mutations cause heritable variations, which are naturally selected, leading to evolution. However, this mutation-led evolution (MLE) is being questioned by an alternative theory called plasticity-led evolution (PLE). PLE suggests that an environmental change induces adaptive phenotypes, which are later genetically accommodated. According to PLE, developmental systems should be able to respond to environmental changes adaptively. However, developmental systems are known to be robust against environmental and mutational perturbations. Thus, we expect a transition from a robust state to a plastic one. To test this hypothesis, we constructed a gene regulatory network (GRN) model that integrates developmental processes, hierarchical regulation, and environmental cues. We then simulated its evolution over different magnitudes of environmental changes. Our findings indicate that this GRN model exhibits PLE under large environmental changes and MLE under small environmental changes. Furthermore, we observed that the GRN model is susceptible to environmental or genetic fluctuations under large environmental changes but is robust under small environmental changes. This indicates a breakdown of robustness due to large environmental changes. Before the breakdown of robustness, the distribution of phenotypes is biased and aligned to the environmental changes, which would facilitate rapid adaptation should a large environmental change occur. These observations suggest that the evolutionary transition from mutation-led to plasticity-led evolution is due to a developmental transition from robust to susceptible regimes over increasing magnitudes of environmental change. Thus, the GRN model can reconcile these conflicting theories of evolution.

Set stabilization of logical control networks: A minimum node control approach

In network systems, control using minimum nodes or pinning control can be effectively used for stabilization problems to cut down the cost of control. In this paper, we investigate the set stabilization problem of logical control networks. In particular, we study the set stabilization problem of probabilistic Boolean networks (PBNs) and probabilistic Boolean control networks (PBCNs) via controlling minimal nodes. Firstly, an algorithm is given to search for the minimum index set of pinning nodes. Then, based on the analysis of its high computational complexity, we present optimized algorithms with lower computational complexity to ascertain the network control using minimum node sets. Moreover, some sufficient and necessary conditions are proposed to ensure the feasibility and effectiveness of the proposed algorithms. Furthermore, a theorem is presented for PBCNs to devise all state-feedback controllers corresponding to the set of pinning nodes. Finally, two models of gene regulatory networks are considered to show the efficacy of obtained results.

RsaL is a self-regulatory switch that controls alternative biosynthesis of two AHL-type quorum sensing signals in Pseudomonas aeruginosa PA1201.

Pseudomonas aeruginosa is a ubiquitous and metabolically versatile microorganism naturally found in soil and water. It is also an opportunistic pathogen in plants, insects, animals, and humans. In response to increasing cell density, P. aeruginosa uses two acyl-homoserine lactone (AHL) quorum-sensing (QS) signals (i.e.,N-3-oxo-dodecanoyl homoserine lactone [3-oxo-C12-HSL] and N-butanoyl-homoserine lactone [C4-HSL]), which regulate the expression of hundreds of genes. However, how the biosynthesis of these two QS signals is coordinated remains unknown. We studied the regulation of these two QS signals in the rhizosphere strain PA1201. PA1201 sequentially produced 3-oxo-C12-HSL and C4-HSL at the early and late growth stages, respectively. The highest 3-oxo-C12-HSL-dependent elastase activity was observed at the early stage, while the highest C4-HSL-dependent rhamnolipid production was observed at the late stage. The atypical regulator RsaL played a pivotal role in coordinating 3-oxo-C12-HSL and C4-HSL biosynthesis and QS-associated virulence. RsaL repressed lasI transcription by binding the -10 and -35 boxes of the lasI promoter. In contrast, RsaL activated rhlI transcription by binding the region encoding the 5'-untranslated region of the rhlI mRNA. Further, RsaL repressed its own expression by binding a nucleotide motif located in the -35 box of the rsaL promoter. Thus, RsaL acts as a molecular switch that coordinates the sequential biosynthesis of AHL QS signals and differential virulence in PA1201. Finally, C4-HSL activation by RsaL was independent of the Las and Pseudomonas quinolone signal (PQS) QS signaling systems. Therefore, we propose a new model of the QS regulatory network in PA1201, in which RsaL represents a superior player acting at the top of the hierarchy.

Chaos in gene regulatory networks: Effects of time delays and interaction structure.

In biological system models, gene expression levels are typically described by regulatory feedback mechanisms. Many studies of gene network models focus on dynamical interactions between components, but often overlook time delays. Here we present an extended model for gene regulatory networks with time delayed negative feedback, which is described by delay differential equations. We analyze nonlinear properties of the model in terms of chaos and compare the conditions with the benchmark homogeneous gene regulatory network model. Chaotic dynamics depend strongly on the inclusion of time delays, but the minimum motifs that show chaos differ when both original and extended models are considered. Our results suggest that, for a particular higher order extension of the gene network, it is possible to observe chaotic dynamics in a two-gene system without adding any self-inhibition. This finding can be explained as a result of modification of the original benchmark model induced by previously unmodeled dynamics. We argue that the inclusion of additional parameters in regulatory gene circuit models substantially enhances the likelihood of observing non-periodic dynamics.

A Spatial Artificial Chemistry Implementation of a Gene Regulatory Network Aimed at Generating Protein Concentration Dynamics.

Gene regulatory networks are networks of interactions in organisms responsible for determining the production levels of proteins and peptides. Mathematical and computational models of gene regulatory networks have been proposed, some of them rather abstract and called artificial regulatory networks. In this contribution, a spatial model for gene regulatory networks is proposed that is biologically more realistic and incorporates an artificial chemistry to realize the interaction between regulatory proteins called the transcription factors and the regulatory sites of simulated genes. The result is a system that is quite robust while able to produce complex dynamics similar to what can be observed in nature. Here an analysis of the impact of the initial states of the system on the produced dynamics is performed, showing that such models are evolvable and can be directed toward producing desired protein dynamics.

A gene regulatory network model that recovers the abaxial-adaxial polarity in Arabidopsis thaliana leaf primordium

Megaphylls, present in the majority of vascular plants, show in many plant lineages an abaxial-adaxial polarity in their dorsoventral axis. This polarity commonly translates into different tissues developing on each side of the leaf blade. This is important because it promotes better photosynthetic efficiency as related to light absorption and gas exchange. Many researchers have studied the molecular bases of the emergence of leaf abaxial-adaxial polarity, showing that it is produced by the interaction and differential expression of particular genes and other molecules. However, until now, it is still unclear if the molecular components documented thus far are sufficient to explain the emergence of leaf polarity. In this work, we integrated the available experimental data to construct a graph of the Gene Regulatory Network (GRN) involved in the formation of abaxial-adaxial polarity in the leaf primordium of Arabidopsis thaliana. This graph consisted of 21 nodes and 47 regulations. We extracted the main components of the graph to obtain a Minimum Network consisting of six genes and 22 possible regulations. Then, we used the Boolean network (BN) formalism to describe the dynamics of this Minimum Network. We identified 1905 distinct BNs that comprised the regulations of the Minimum Network and exclusively generated the two attractors representing the abaxial and adaxial cell types. This highlights the fact that most graphs, including our network, can describe experimentally observed behaviors with many BN dynamics. By performing mutant simulations and robustness analysis, we found that two of the 1905 BNs better reproduce experimentally available information. To produce the expected attractors, both BNs predict the same missing regulations, which we propose should be experimentally analyzed to confirm their existence. Interestingly, these two BNs have low robustness to perturbations compared with previously analyzed GRNs. This was an unexpected result since abaxial-adaxial polarity is a robust biological trait, which suggests more components or regulations of the network are missing.

A Kinetic Finite Volume Discretization of the Multidimensional PIDE Model for Gene Regulatory Networks

In this paper, a finite volume discretization scheme for partial integro-differential equations (PIDEs) describing the temporal evolution of protein distribution in gene regulatory networks is proposed. It is shown that the obtained set of ODEs can be formally represented as a compartmental kinetic system with a strongly connected reaction graph. This allows the application of the theory of nonnegative and compartmental systems for the qualitative analysis of the approximating dynamics. In this framework, it is straightforward to show the existence, uniqueness and stability of equilibria. Moreover, the computation of the stationary probability distribution can be traced back to the solution of linear equations. The discretization scheme is presented for one and multiple dimensional models separately. Illustrative computational examples show the precision of the approach, and good agreement with previous results in the literature.