In this thesis, we examine representation of positive integers by certain definite quaternary quadratic forms Q over ZZ and ZZ [(1+ sqrt{5})/2] by studying the theta series of the forms, which are (Hilbert) modular forms with weight 2 and of level and character determined by Q. Specifically, we write the theta series for the following three quadratic forms as the sum of an explicit Eisenstein series and linear combination of normalized Hecke eigencusp forms: Q_1(x) = x_1^2+x_2^2+x_3^2+x_4^2 over ZZ (which has long been studied, and which we provide for the sake of introduction), Q_2(x) = x_1^2+x_2^2+x_3^2+7x_4^2 over ZZ (which seems not to have appeared this explicitly previously in the literature), and of Q_1 over ZZ[(1+ sqrt{5})/2]. We also provide an explicit formula for the Eisenstein series Fourier coefficients appearing in the theta series associated to Q_3(x) = x_1^2+x_2^2+3x_3^2+3(3+sqrt{5})x_4^2 over ZZ[(1+ sqrt{5})/2]. Beyond representation results, we develop and implement an algorithm which counts the number of representations of totally positive integers by a quaternary definite integral quadratic form Q defined over the ring of integers of any real quadratic number field. We also implement an algorithm of Dembele which returns a Hecke eigenbasis for spaces of Hilbert modular cusp forms of parallel weight two and trivial character via an explicit version of the Eichler-Shimizu-Jacquet-Langlands correspondence.A primary tool for studying the theta series is Siegel's product formula to compute the Eisenstein components of the Fourier coefficients. To understand the cuspidal component of the theta series, we first find a Hecke eigenbasis for the corresponding modular cusp space and write the cuspidal component of the theta series as a linear combination of these basis elements. Then, using known asymptotics on the Eisenstein coefficients and bounds by Deligne on the cusp form coefficients we prove that all sufficiently large locally represented integers are represented. In particular for the form Q_2, we show that all locally represented m >15825811 are represented; we also present for this form upper and lower bounds on r_Q(m) explicitly obtained from the theta series decomposition.