A machine learning-based constitutive model is proposed to capture the strain rate-dependent, nonlinear stress-strain response of soft materials. This model employs two Gaussian process regression (GPR) models to discover the complex constitutive relationships between stress, strain, and strain rate in the material. Specifically, the first GPR captures the elastic stress component by mapping invariants of the right Cauchy-Green deformation tensor, to a set of coefficients that define the elastic first Piola-Kirchhoff stress tensor as a linear combination of an irreducible integrity basis. The second GPR captures the viscous overstress by mapping the invariants of the right Cauchy-Green deformation tensor and its rate of change with respect to time, to the corresponding set of coefficients of the viscous Piola-Kirchhoff overstress tensor through a similar procedure. It is shown that this type of model construction automatically ensures a number of physical constraints on material behavior: frame-indifference, symmetry, balance of angular momentum, stress-free initial state, principles of determinism and local action, and second law of thermodynamics for the elastic stress component. Additionally, it is shown to efficiently generalize beyond the training region and perform efficiently with limited data. To further impose thermodynamic consistency of the viscous overstress component, the log-marginal likelihood function of the second GPR is optimized under an inequality constraint that requires that the tensor contraction of the viscous overstress with the rate of the right Cauchy-Green deformation tensor is strictly non-negative. The two GPR models are then trained using artificially generated stress-strain-strain rate datasets from available hyperelastic and visco-hyperelastic models, respectively. From preliminary results of the elastic stress GPR model, it is seen that the data-driven constitutive model offers both accurate fitting and a good predictive accuracy.