A common task in experimental sciences is to fit mathematical models to real-world measurements to improve understanding of natural phenomenon (reverse-engineering or inverse modeling). When complex dynamical systems are considered, such as partial differential equations, this task may become challenging and ill-posed. In this work, a linear parabolic equation is considered where the objective is to estimate both the differential operator coefficients and the source term at once. The Bayesian methodology for inverse problems provides a form of regularization while quantifying uncertainty as the solution is a probability measure taking in account data. This posterior distribution, which is non-Gaussian and infinite dimensional, is then summarized through a mode and sampled using a state-of-the-art Markov-Chain Monte-Carlo algorithm based on a geometric approach. After a rigorous analysis, this methodology is applied on a dataset of the post-transcriptional regulation of Kni gap gene in the early development of Drosophila Melanogaster where mRNA concentration and both diffusion and depletion rates are inferred from noisy measurement of the protein concentration