Computational systems biology is an emerging field built upon advances in biological experimentation, computational tools, and joint interdisciplinary efforts from physical sciences to engineering. In this thesis, I solve two important problems in the developmental genetics of Drosophila melanogaster by means of two different computational systems approaches, one stochastic and the other deterministic. In the first problem I use exact simulations of a spatially extended Master Equation to analyze intrinsic fluctuations in the gradient of the Bicoid morphogen. These simulation results were in turn compared to quantitative experimental measurements of Bicoid levels performed in fixed tissue by immunostaining. We selected Bicoid for this study because different concentrations of Bicoid can trigger different developmental pathways and determine the cell fate. These properties define a morphogen gradient, and Bicoid is the best characterized such gradient at the molecular level. At the theoretical level, the Bicoid gradient is amenable to investigation by stochastic models because its formation can be described by the processes of diffusion and decay, with synthesis described by a boundary condition. I show that the intrinsic noise of Bicoid gradient is Poisson distributed. I demonstrate how experimental noise can be identified in the logarithm domain from single embryo analysis, and then separated from intrinsic noise in the normalized variance domain of an ensemble statistical analysis. I show how measurement sensitivity affects our observations, and how small amounts of rescaling noise can perturb the noise strength (Fano factor) observed. I demonstrate that the biological noise level in data can serve as a physical constraint for restricting the model's parameter space, and for predicting the Bicoid molecular number and variation range. I exhibit the predicted molecular number gradient together with measurement effects, and make a comparison between conditions of higher and lower variance respectively. In the second problem I use the gene circuit method, which uses coarse-grained chemical kinetic equations fitted by a data-driven optimization approach, to study the pattern formation regulatory control of pair-rule genes. The mathematical model, based on mechanisms of protein synthesis, decay and diffusion, captures the essential control of the size, location and timing of the pair-rule stripe formation process. The gene circuits were optimized in modular manner, such that two classes of circuits were found by systematically and combinatorially reducing optimization constraints. Three major circuits were selected for detailed dynamical regulatory analysis. These circuits include 5 cross-regulating pair rule genes (eve, h, run, ftz, odd) and 7 external input gap and maternal genes (bcd, cad, hb, Kr, gt, kni, tll). The circuits model a period of time from cycle 13 to gastrulation in a spatial region covering 44 nuclei from 35% to 78% EL, extending from eve stripe 2 to the stripe 7 anterior border. The biological conclusions were drawn in 6 levels based on a consensus among circuits and literature. The regulatory phasing analysis examines a simple phasing rule based on stripe formation and shifting constraints. The major predictions include gap gene activation on odd, with exceptionally strong consensus among all circuits. The complete lack of understanding and literature reference for gap gene regulations of odd illustrates the significance and advantage of the computational system approach. Other systems conclusions include that all gap gene regulation on eve is repressive, while all gap gene regulation on odd is activating. h has the least gap gene inputs in the circuits, which is strong contradiction to literature.