A soft tissue’s macroscopic behavior is largely determined by its microstructural components (often a collagen fiber network surrounded by a non-fibrillar matrix (NFM)). In the present study, a coupled fiber-matrix model was developed to quantify fully the internal stress field within such a tissue and to explore interactions between the collagen fiber network and non-fibrillar matrix (NFM). Voronoi tessellations (representing collagen networks) were embedded in a continuous three-dimensional NFM. Fibers were represented as one-dimensional nonlinear springs and the NFM, meshed via tetrahedra, was modeled as a compressible neo-Hookean solid. Multidimensional finite element modeling was employed to couple the two tissue components, and uniaxial tension was applied to the composite representative volume element (RVE). In terms of the overall RVE response (average stress, fiber orientation, Poisson’s ratio), the coupled fiber-matrix model yielded results consistent with those obtained using a previously developed parallel model based upon superposition. The detailed stress field in the composite RVE demonstrated the high degree of in homogeneity in NFM mechanics, which cannot be addressed by a parallel model. Distributions of maximum/minimum principal stresses in the NFM showed a transition from fiber-dominated to matrix-dominated behavior as the matrix shear modulus increased. The matrix-dominated behavior also included a shift in the fiber kinematics toward the affine limit. We conclude that if only gross averaged parameters are of interest, parallel-type models are suitable. If, however, one is concerned with phenomena, such as individual cell-fiber interactions or tissue failure, that could be altered by local variations in the stress field, then the detailed model is necessary in spite of its higher computational cost.