Representative single anisotropic matrix block 2D Green’s function models for depletion through fully-penetrating, vertical fractures through different numbers of fracture faces are constructed that analytically capture both fracture and block depletion with fracture-matrix mass transfer. The 1D Green’s function for a fracture system is likewise solved in terms of the time evolution of average fracture pressure. While transient average pressure values are not inherently measurable, they are transformed into cumulative production or instantaneous flowrate values, thus producing new rate decline model functional forms. Primary variables in assembling the interacting systems model are the volume ratio, Vf /Vm, permeability ratio, kf /kx, and geometry, (a/b)(ky/kx), with the last term accounting for both block shape and permeability anisotropy. We construct interacting systems models in terms of various ratios of Vf /Vm, and kf /kx for three fracture architecture prototypes: representative matrix blocks depleted by 4, 2, or 1 contacting fractures. The single matrix block models can be migrated to ones for heterogeneous systems using superposition and matrix block distributions, as demonstrated with a binary distribution of block sizes with variable fractions. Analytic solutions for rate decline problems can be used to understand the production signatures of naturally fractured reservoirs and interpretation of fracture volume fraction, permeability ratio, average matrix block size, and measures of heterogeneity.