51
|
Pozo C, Marín-Sanguino A, Alves R, Guillén-Gosálbez G, Jiménez L, Sorribas A. Steady-state global optimization of metabolic non-linear dynamic models through recasting into power-law canonical models. BMC SYSTEMS BIOLOGY 2011; 5:137. [PMID: 21867520 PMCID: PMC3201032 DOI: 10.1186/1752-0509-5-137] [Citation(s) in RCA: 20] [Impact Index Per Article: 1.4] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Subscribe] [Scholar Register] [Received: 04/05/2011] [Accepted: 08/25/2011] [Indexed: 01/18/2023]
Abstract
Background Design of newly engineered microbial strains for biotechnological purposes would greatly benefit from the development of realistic mathematical models for the processes to be optimized. Such models can then be analyzed and, with the development and application of appropriate optimization techniques, one could identify the modifications that need to be made to the organism in order to achieve the desired biotechnological goal. As appropriate models to perform such an analysis are necessarily non-linear and typically non-convex, finding their global optimum is a challenging task. Canonical modeling techniques, such as Generalized Mass Action (GMA) models based on the power-law formalism, offer a possible solution to this problem because they have a mathematical structure that enables the development of specific algorithms for global optimization. Results Based on the GMA canonical representation, we have developed in previous works a highly efficient optimization algorithm and a set of related strategies for understanding the evolution of adaptive responses in cellular metabolism. Here, we explore the possibility of recasting kinetic non-linear models into an equivalent GMA model, so that global optimization on the recast GMA model can be performed. With this technique, optimization is greatly facilitated and the results are transposable to the original non-linear problem. This procedure is straightforward for a particular class of non-linear models known as Saturable and Cooperative (SC) models that extend the power-law formalism to deal with saturation and cooperativity. Conclusions Our results show that recasting non-linear kinetic models into GMA models is indeed an appropriate strategy that helps overcoming some of the numerical difficulties that arise during the global optimization task.
Collapse
Affiliation(s)
- Carlos Pozo
- Departament de Ciències Mèdiques Bàsiques, Institut de Recerca Biomèdica de Lleida (IRBLLEIDA), Universitat de Lleida, Spain
| | | | | | | | | | | |
Collapse
|
52
|
Schmidt MD, Vallabhajosyula RR, Jenkins JW, Hood JE, Soni AS, Wikswo JP, Lipson H. Automated refinement and inference of analytical models for metabolic networks. Phys Biol 2011; 8:055011. [PMID: 21832805 DOI: 10.1088/1478-3975/8/5/055011] [Citation(s) in RCA: 45] [Impact Index Per Article: 3.2] [Reference Citation Analysis] [Abstract] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 01/25/2023]
Abstract
The reverse engineering of metabolic networks from experimental data is traditionally a labor-intensive task requiring a priori systems knowledge. Using a proven model as a test system, we demonstrate an automated method to simplify this process by modifying an existing or related model--suggesting nonlinear terms and structural modifications--or even constructing a new model that agrees with the system's time series observations. In certain cases, this method can identify the full dynamical model from scratch without prior knowledge or structural assumptions. The algorithm selects between multiple candidate models by designing experiments to make their predictions disagree. We performed computational experiments to analyze a nonlinear seven-dimensional model of yeast glycolytic oscillations. This approach corrected mistakes reliably in both approximated and overspecified models. The method performed well to high levels of noise for most states, could identify the correct model de novo, and make better predictions than ordinary parametric regression and neural network models. We identified an invariant quantity in the model, which accurately derived kinetics and the numerical sensitivity coefficients of the system. Finally, we compared the system to dynamic flux estimation and discussed the scaling and application of this methodology to automated experiment design and control in biological systems in real time.
Collapse
Affiliation(s)
- Michael D Schmidt
- Cornell Computational Systems Laboratory, Cornell University, Ithaca, NY, USA
| | | | | | | | | | | | | |
Collapse
|
53
|
Costa R, Rocha I, Ferreira E, Machado D. Critical perspective on the consequences of the limited availability of kinetic data in metabolic dynamic modelling. IET Syst Biol 2011; 5:157-63. [DOI: 10.1049/iet-syb.2009.0058] [Citation(s) in RCA: 18] [Impact Index Per Article: 1.3] [Reference Citation Analysis] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 12/18/2022] Open
|
54
|
Zhan C, Yeung LF. Parameter estimation in systems biology models using spline approximation. BMC SYSTEMS BIOLOGY 2011; 5:14. [PMID: 21255466 PMCID: PMC3750107 DOI: 10.1186/1752-0509-5-14] [Citation(s) in RCA: 53] [Impact Index Per Article: 3.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Subscribe] [Scholar Register] [Received: 03/18/2010] [Accepted: 01/24/2011] [Indexed: 11/24/2022]
Abstract
Background Mathematical models for revealing the dynamics and interactions properties of biological systems play an important role in computational systems biology. The inference of model parameter values from time-course data can be considered as a "reverse engineering" process and is still one of the most challenging tasks. Many parameter estimation methods have been developed but none of these methods is effective for all cases and can overwhelm all other approaches. Instead, various methods have their advantages and disadvantages. It is worth to develop parameter estimation methods which are robust against noise, efficient in computation and flexible enough to meet different constraints. Results Two parameter estimation methods of combining spline theory with Linear Programming (LP) and Nonlinear Programming (NLP) are developed. These methods remove the need for ODE solvers during the identification process. Our analysis shows that the augmented cost function surfaces used in the two proposed methods are smoother; which can ease the optima searching process and hence enhance the robustness and speed of the search algorithm. Moreover, the cores of our algorithms are LP and NLP based, which are flexible and consequently additional constraints can be embedded/removed easily. Eight system biology models are used for testing the proposed approaches. Our results confirm that the proposed methods are both efficient and robust. Conclusions The proposed approaches have general application to identify unknown parameter values of a wide range of systems biology models.
Collapse
Affiliation(s)
- Choujun Zhan
- Department of Electronic Engineering, City University of Hong Kong, PR China.
| | | |
Collapse
|
55
|
Garcia J, Sims KJ, Schwacke JH, Del Poeta M. Biochemical systems analysis of signaling pathways to understand fungal pathogenicity. Methods Mol Biol 2011; 734:173-200. [PMID: 21468990 PMCID: PMC5155339 DOI: 10.1007/978-1-61779-086-7_9] [Citation(s) in RCA: 6] [Impact Index Per Article: 0.4] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Indexed: 02/12/2023]
Abstract
Over the past decade, researchers have recognized the need to study biological systems as integrated systems. While the reductionist approaches of the past century have made remarkable advances of our understanding of life, the next phase of understanding comes from systems-level investigations. Additionally, biology has become a data-intensive field of research. The introduction of high throughput sequencing, microarrays, high throughput proteomics, metabolomics, and now lipidomics are producing significantly more data than can be interpreted using existing methods. The field of systems biology brings together methods from computer science, modeling, statistics, engineering, and biology to explore the volumes of data now being produced and to develop mathematical representations of metabolic, signaling, and gene regulatory systems. Advances in these methods are allowing biologists to develop new insights into the complexities of life, to predict cellular responses and treatment outcomes, and to effectively plan experiments that extend our understanding. In this chapter, we are providing the basic steps of developing and analyzing a small S-system model of a biochemical pathway related to sphingolipid metabolism in the regulation of virulence of the human fungal microbial pathogen Cryptococcus neoformans (Cn).
Collapse
Affiliation(s)
- Jacqueline Garcia
- Department of Biochemistry and Molecular Biology, Medical University of South Carolina, Charleston, South Carolina, USA
| | - Kellie J Sims
- Department of Biochemistry and Molecular Biology, Medical University of South Carolina, Charleston, South Carolina, USA
| | - John H Schwacke
- Department of Biochemistry and Molecular Biology, Medical University of South Carolina, Charleston, South Carolina, USA
| | - Maurizio Del Poeta
- Department of Biochemistry and Molecular Biology, Medical University of South Carolina, Charleston, South Carolina, USA
- Department of Microbiology and Immunology, Medical University of South Carolina, Charleston, South Carolina, USA
- Department of Craniofacial Biology, Medical University of South Carolina, Charleston, South Carolina, USA
- Division of Infectious Disease, Medical University of South Carolina, Charleston, South Carolina, USA
| |
Collapse
|
56
|
Fonseca LL, Sánchez C, Santos H, Voit EO. Complex coordination of multi-scale cellular responses to environmental stress. MOLECULAR BIOSYSTEMS 2010; 7:731-41. [PMID: 21088798 DOI: 10.1039/c0mb00102c] [Citation(s) in RCA: 17] [Impact Index Per Article: 1.1] [Reference Citation Analysis] [Abstract] [Track Full Text] [Subscribe] [Scholar Register] [Indexed: 11/21/2022]
Abstract
Cells and organisms are regularly exposed to a variety of stresses, and effective responses are a matter of survival. The article describes a multi-scale experimental and dynamical modeling analysis that clearly indicates concerted stress control in different temporal and organizational domains, and a strong synergy between the dynamics of genes, proteins and metabolites. Specifically, we show with in vivo NMR measurements of metabolic profiles that baker's yeast responds to a paradigmatic stress, heat, at three organizational levels and in two time regimes. At the metabolic level, an almost immediate response is mounted. However, this response is a "quick fix" in comparison to a much more effective response that had been pre-organized in earlier periods of heat stress and is an order of magnitude stronger. Equipped with the metabolic profile data, our modeling efforts resulted in a crisp, quantitative separation of response actions at the levels of metabolic control and gene regulation. They also led to predictions of necessary changes in protein levels and clearly demonstrated that formerly observed temperature profiles of key enzyme activities are not sufficient to explain the accumulation of trehalose as an immediate response to sudden heat stress.
Collapse
Affiliation(s)
- Luís L Fonseca
- Instituto de Tecnologia Química e Biológica, Universidade Nova de Lisboa, 2780-156 Oeiras, Portugal
| | | | | | | |
Collapse
|
57
|
Guillén-Gosálbez G, Sorribas A. Identifying quantitative operation principles in metabolic pathways: a systematic method for searching feasible enzyme activity patterns leading to cellular adaptive responses. BMC Bioinformatics 2009; 10:386. [PMID: 19930714 PMCID: PMC2799421 DOI: 10.1186/1471-2105-10-386] [Citation(s) in RCA: 25] [Impact Index Per Article: 1.6] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 02/25/2009] [Accepted: 11/24/2009] [Indexed: 01/13/2023] Open
Abstract
BACKGROUND Optimization methods allow designing changes in a system so that specific goals are attained. These techniques are fundamental for metabolic engineering. However, they are not directly applicable for investigating the evolution of metabolic adaptation to environmental changes. Although biological systems have evolved by natural selection and result in well-adapted systems, we can hardly expect that actual metabolic processes are at the theoretical optimum that could result from an optimization analysis. More likely, natural systems are to be found in a feasible region compatible with global physiological requirements. RESULTS We first present a new method for globally optimizing nonlinear models of metabolic pathways that are based on the Generalized Mass Action (GMA) representation. The optimization task is posed as a nonconvex nonlinear programming (NLP) problem that is solved by an outer-approximation algorithm. This method relies on solving iteratively reduced NLP slave subproblems and mixed-integer linear programming (MILP) master problems that provide valid upper and lower bounds, respectively, on the global solution to the original NLP. The capabilities of this method are illustrated through its application to the anaerobic fermentation pathway in Saccharomyces cerevisiae. We next introduce a method to identify the feasibility parametric regions that allow a system to meet a set of physiological constraints that can be represented in mathematical terms through algebraic equations. This technique is based on applying the outer-approximation based algorithm iteratively over a reduced search space in order to identify regions that contain feasible solutions to the problem and discard others in which no feasible solution exists. As an example, we characterize the feasible enzyme activity changes that are compatible with an appropriate adaptive response of yeast Saccharomyces cerevisiae to heat shock CONCLUSION Our results show the utility of the suggested approach for investigating the evolution of adaptive responses to environmental changes. The proposed method can be used in other important applications such as the evaluation of parameter changes that are compatible with health and disease states.
Collapse
Affiliation(s)
- Gonzalo Guillén-Gosálbez
- Departament de Ciències Mèdiques Bàsiques, Institut de Recerca Biomèdica de Lleida, Universitat de Lleida, Montserrat Roig 2, 25008-Lleida, Spain.
| | | |
Collapse
|
58
|
Yin W, Jo H, Voit EO. Systems analysis of the role of bone morphogenic protein 4 in endothelial inflammation. Ann Biomed Eng 2009; 38:291-307. [PMID: 19851868 DOI: 10.1007/s10439-009-9822-y] [Citation(s) in RCA: 9] [Impact Index Per Article: 0.6] [Reference Citation Analysis] [Abstract] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 04/08/2009] [Accepted: 10/13/2009] [Indexed: 02/07/2023]
Abstract
Shear stress is an important factor in the onset and progression of atherosclerosis. High and unidirectional laminar stress is seen as protective, while low and oscillatory shear stress is considered pro-inflammatory and pro-atherogenic. The mechanosensitive response of endothelial cells is governed by a complex system of genes, proteins, and signals that operate at distinctly different time scales. We propose a dynamic mathematical model that quantitatively describes this mechanosensing system and permits novel insights into its functioning. The model, the first of its kind, is constructed within the guidelines of Biochemical Systems Theory and accounts for different time scales by means of approximated delays. Parameter values are obtained directly from biochemical observations in an ad hoc fashion. The model reflects most documented observations well and leads to a number of predictions and novel hypotheses. In particular, it demonstrates the crucial role of Bone Morphogenic Protein 4 and p47(phox)-dependent NADPH oxidases in endothelial inflammation.
Collapse
Affiliation(s)
- Weiwei Yin
- Integrative BioSystems Institute, Georgia Institute of Technology, 313 Ferst Drive, Atlanta, GA 30332, USA.
| | | | | |
Collapse
|
59
|
Calibration of dynamic models of biological systems with KInfer. EUROPEAN BIOPHYSICS JOURNAL: EBJ 2009; 39:1019-39. [PMID: 19669750 DOI: 10.1007/s00249-009-0520-3] [Citation(s) in RCA: 8] [Impact Index Per Article: 0.5] [Reference Citation Analysis] [Abstract] [Track Full Text] [Subscribe] [Scholar Register] [Received: 02/16/2009] [Revised: 07/02/2009] [Accepted: 07/05/2009] [Indexed: 02/03/2023]
Abstract
Methods for parameter estimation that are robust to experimental uncertainties and to stochastic and biological noise and that require a minimum of a priori input knowledge are of key importance in computational systems biology. The new method presented in this paper aims to ensure an inference model that deduces the rate constants of a system of biochemical reactions from experimentally measured time courses of reactants. This new method was applied to some challenging parameter estimation problems of nonlinear dynamic biological systems and was tested both on synthetic and real data. The synthetic case studies are the 12-state model of the SERCA pump and a model of a genetic network containing feedback loops of interaction between regulator and effector genes. The real case studies consist of a model of the reaction between the inhibitor kappaB kinase enzyme and its substrate in the signal transduction pathway of NF-kappaB, and a stiff model of the fermentation pathway of Lactococcus lactis.
Collapse
|
60
|
Chou IC, Voit EO. Recent developments in parameter estimation and structure identification of biochemical and genomic systems. Math Biosci 2009; 219:57-83. [PMID: 19327372 PMCID: PMC2693292 DOI: 10.1016/j.mbs.2009.03.002] [Citation(s) in RCA: 298] [Impact Index Per Article: 18.6] [Reference Citation Analysis] [Abstract] [Key Words] [MESH Headings] [Grants] [Track Full Text] [Journal Information] [Subscribe] [Scholar Register] [Received: 10/02/2008] [Revised: 03/06/2009] [Accepted: 03/15/2009] [Indexed: 01/16/2023]
Abstract
The organization, regulation and dynamical responses of biological systems are in many cases too complex to allow intuitive predictions and require the support of mathematical modeling for quantitative assessments and a reliable understanding of system functioning. All steps of constructing mathematical models for biological systems are challenging, but arguably the most difficult task among them is the estimation of model parameters and the identification of the structure and regulation of the underlying biological networks. Recent advancements in modern high-throughput techniques have been allowing the generation of time series data that characterize the dynamics of genomic, proteomic, metabolic, and physiological responses and enable us, at least in principle, to tackle estimation and identification tasks using 'top-down' or 'inverse' approaches. While the rewards of a successful inverse estimation or identification are great, the process of extracting structural and regulatory information is technically difficult. The challenges can generally be categorized into four areas, namely, issues related to the data, the model, the mathematical structure of the system, and the optimization and support algorithms. Many recent articles have addressed inverse problems within the modeling framework of Biochemical Systems Theory (BST). BST was chosen for these tasks because of its unique structural flexibility and the fact that the structure and regulation of a biological system are mapped essentially one-to-one onto the parameters of the describing model. The proposed methods mainly focused on various optimization algorithms, but also on support techniques, including methods for circumventing the time consuming numerical integration of systems of differential equations, smoothing overly noisy data, estimating slopes of time series, reducing the complexity of the inference task, and constraining the parameter search space. Other methods targeted issues of data preprocessing, detection and amelioration of model redundancy, and model-free or model-based structure identification. The total number of proposed methods and their applications has by now exceeded one hundred, which makes it difficult for the newcomer, as well as the expert, to gain a comprehensive overview of available algorithmic options and limitations. To facilitate the entry into the field of inverse modeling within BST and related modeling areas, the article presented here reviews the field and proposes an operational 'work-flow' that guides the user through the estimation process, identifies possibly problematic steps, and suggests corresponding solutions based on the specific characteristics of the various available algorithms. The article concludes with a discussion of the present state of the art and with a description of open questions.
Collapse
Affiliation(s)
- I-Chun Chou
- Integrative BioSystems Institute and The Wallace H. Coulter Department of Biomedical Engineering, Georgia Institute of Technology and Emory University, 313 Ferst Drive, Atlanta, GA 30332, USA.
| | | |
Collapse
|
61
|
Ko CL, Voit EO, Wang FS. Estimating parameters for generalized mass action models with connectivity information. BMC Bioinformatics 2009; 10:140. [PMID: 19432964 PMCID: PMC2694188 DOI: 10.1186/1471-2105-10-140] [Citation(s) in RCA: 12] [Impact Index Per Article: 0.8] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Journal Information] [Subscribe] [Scholar Register] [Received: 02/03/2009] [Accepted: 05/11/2009] [Indexed: 11/10/2022] Open
Abstract
BACKGROUND Determining the parameters of a mathematical model from quantitative measurements is the main bottleneck of modelling biological systems. Parameter values can be estimated from steady-state data or from dynamic data. The nature of suitable data for these two types of estimation is rather different. For instance, estimations of parameter values in pathway models, such as kinetic orders, rate constants, flux control coefficients or elasticities, from steady-state data are generally based on experiments that measure how a biochemical system responds to small perturbations around the steady state. In contrast, parameter estimation from dynamic data requires time series measurements for all dependent variables. Almost no literature has so far discussed the combined use of both steady-state and transient data for estimating parameter values of biochemical systems. RESULTS In this study we introduce a constrained optimization method for estimating parameter values of biochemical pathway models using steady-state information and transient measurements. The constraints are derived from the flux connectivity relationships of the system at the steady state. Two case studies demonstrate the estimation results with and without flux connectivity constraints. The unconstrained optimal estimates from dynamic data may fit the experiments well, but they do not necessarily maintain the connectivity relationships. As a consequence, individual fluxes may be misrepresented, which may cause problems in later extrapolations. By contrast, the constrained estimation accounting for flux connectivity information reduces this misrepresentation and thereby yields improved model parameters. CONCLUSION The method combines transient metabolic profiles and steady-state information and leads to the formulation of an inverse parameter estimation task as a constrained optimization problem. Parameter estimation and model selection are simultaneously carried out on the constrained optimization problem and yield realistic model parameters that are more likely to hold up in extrapolations with the model.
Collapse
Affiliation(s)
- Chih-Lung Ko
- Department of Chemical Engineering, National Chung Cheng University, Chiayi, Taiwan.
| | | | | |
Collapse
|
62
|
Vilela M, Vinga S, Maia MAGM, Voit EO, Almeida JS. Identification of neutral biochemical network models from time series data. BMC SYSTEMS BIOLOGY 2009. [PMID: 19416537 DOI: 10.1186/1752‐0509‐3‐47] [Citation(s) in RCA: 1] [Impact Index Per Article: 0.1] [Reference Citation Analysis] [Abstract] [Subscribe] [Scholar Register] [Indexed: 11/10/2022]
Abstract
BACKGROUND The major difficulty in modeling biological systems from multivariate time series is the identification of parameter sets that endow a model with dynamical behaviors sufficiently similar to the experimental data. Directly related to this parameter estimation issue is the task of identifying the structure and regulation of ill-characterized systems. Both tasks are simplified if the mathematical model is canonical, i.e., if it is constructed according to strict guidelines. RESULTS In this report, we propose a method for the identification of admissible parameter sets of canonical S-systems from biological time series. The method is based on a Monte Carlo process that is combined with an improved version of our previous parameter optimization algorithm. The method maps the parameter space into the network space, which characterizes the connectivity among components, by creating an ensemble of decoupled S-system models that imitate the dynamical behavior of the time series with sufficient accuracy. The concept of sloppiness is revisited in the context of these S-system models with an exploration not only of different parameter sets that produce similar dynamical behaviors but also different network topologies that yield dynamical similarity. CONCLUSION The proposed parameter estimation methodology was applied to actual time series data from the glycolytic pathway of the bacterium Lactococcus lactis and led to ensembles of models with different network topologies. In parallel, the parameter optimization algorithm was applied to the same dynamical data upon imposing a pre-specified network topology derived from prior biological knowledge, and the results from both strategies were compared. The results suggest that the proposed method may serve as a powerful exploration tool for testing hypotheses and the design of new experiments.
Collapse
Affiliation(s)
- Marco Vilela
- Instituto de Tecnologia Química e Biológica, Universidade Nova de Lisboa, Apartado, Oeiras, Portugal.
| | | | | | | | | |
Collapse
|
63
|
Vilela M, Vinga S, Maia MAGM, Voit EO, Almeida JS. Identification of neutral biochemical network models from time series data. BMC SYSTEMS BIOLOGY 2009; 3:47. [PMID: 19416537 PMCID: PMC2694766 DOI: 10.1186/1752-0509-3-47] [Citation(s) in RCA: 36] [Impact Index Per Article: 2.3] [Reference Citation Analysis] [Abstract] [Track Full Text] [Download PDF] [Figures] [Subscribe] [Scholar Register] [Received: 01/18/2009] [Accepted: 05/05/2009] [Indexed: 12/02/2022]
Abstract
Background The major difficulty in modeling biological systems from multivariate time series is the identification of parameter sets that endow a model with dynamical behaviors sufficiently similar to the experimental data. Directly related to this parameter estimation issue is the task of identifying the structure and regulation of ill-characterized systems. Both tasks are simplified if the mathematical model is canonical, i.e., if it is constructed according to strict guidelines. Results In this report, we propose a method for the identification of admissible parameter sets of canonical S-systems from biological time series. The method is based on a Monte Carlo process that is combined with an improved version of our previous parameter optimization algorithm. The method maps the parameter space into the network space, which characterizes the connectivity among components, by creating an ensemble of decoupled S-system models that imitate the dynamical behavior of the time series with sufficient accuracy. The concept of sloppiness is revisited in the context of these S-system models with an exploration not only of different parameter sets that produce similar dynamical behaviors but also different network topologies that yield dynamical similarity. Conclusion The proposed parameter estimation methodology was applied to actual time series data from the glycolytic pathway of the bacterium Lactococcus lactis and led to ensembles of models with different network topologies. In parallel, the parameter optimization algorithm was applied to the same dynamical data upon imposing a pre-specified network topology derived from prior biological knowledge, and the results from both strategies were compared. The results suggest that the proposed method may serve as a powerful exploration tool for testing hypotheses and the design of new experiments.
Collapse
Affiliation(s)
- Marco Vilela
- Instituto de Tecnologia Química e Biológica, Universidade Nova de Lisboa, Apartado, Oeiras, Portugal.
| | | | | | | | | |
Collapse
|