Classical rigid musculoskeletal lumbar spine models cannot provide spinal load within intervertebral discs (IVDs) due to modeling assumptions. The present study develops an approach for predicting IVD stress in annulus fibrosus (AF) and nucleus pulposus (NP) regions. A musculoskeletal model coupled with elastic elements placed between vertebrae is developed. These elements are distributed within AF and NP regions of each IVD. Each elastic element representing a standard nonlinear solid model allows vertebral strain-like displacement and IVD stress during flexion motion to be computed. Maximal stresses at the L4-5 level are 1.25 and 1.4 MPa for AF and NP regions, respectively, during 40° flexion motion. Maximal stresses at the L2-3, L1-2, and L3-4 levels are 1.93 and 2.2 MPa, 1.12 and 1.2 MPa, and 2.68 MPa for AF and NP regions, respectively. Vertebral displacements, muscle forces, and predicted stresses are in agreement with data reported in the literature. This study opens new perspectives for further efficient musculoskeletal spine models for predicting spinal loads, leading to improved and optimized treatment, recovery, and follow-up recommendations for spinal disorders.