In recent years, additive manufacturing has shown great potential for the design of geometrically complex, biocompatible, and biodegradable scaffold-assisted implants in bone tissue engineering. Mechanoregulatory models have shown that the differentiation of the cells on the scaffold is impacted by quantities such as the surface shear stress which appears when the scaffold is mechanically loaded, and the wall shear stress that emerges due to the interaction of the scaffold with the blood flow around it. However, the performance of the scaffold depends on various parameters such as material composition, geometry, loading, and flow conditions. In this work, to gain an understanding of the parameters' impact on the behavior of the scaffold, a four-layer/strut orthogonal scaffold with isometric pores made of polylactic acid and reinforced by 5% stainless steel particles, is numerically studied. The high-fidelity structural analysis is performed with the Boundary Element Method while for the fluid simulations, the element-based Finite Volume Method is employed. The effect of the variability in the material properties and the boundary conditions of the computational models is explored by utilizing probabilistic deep learning-based reduced order models (ROMs) for the structural and fluid problem. Regarding the training of the ROMs, it is shown that a low number of full-scale simulations is required to compute the high-fidelity data. To implement the ROMs, proper orthogonal decomposition is used to compute a low dimensional basis, and a Bayesian neural network to map the simulation parameters to the reduced data as well as to capture the epistemic uncertainties in the regression task. The predicted results demonstrate the ability of the ROMs to provide numerical results with good accuracy for the parameter range of interest.