To economically develop shale gas reservoirs, hydraulic fracturing treatments must be applied to artificially introduce pathways for gas to flow out. After hydraulic fracturing, both the induced hydraulic fractures and activated natural fractures provide high permeable channels for gas transport. However, the long-term production performance of shale gas reservoirs is determined by the physical properties of matrix instead of the fractures. In this work, a coupled hydro-mechanical model is developed to investigate gas flow in heterogeneous shale matrix at the microscale. This model considers the mechanical interaction and gas mass transfer between organic matter (OM) and inorganic matter (iOM) during gas depletion. The interesting results show that due to matrix heterogeneity, i.e. co-existence of OM and iOM, the gas flow and mechanical deformation behaviors in shale exhibit regional differences. Specifically, during reservoir depletion, pore pressure, Knudsen number and volumetric strain in iOM are smaller than that in OM while the intrinsic and apparent permeability of OM are smaller than that of iOM. Thus, matrix heterogeneity should be taken into account when characterizing shale gas flow. In addition, sensitivity analysis is conducted to study the influence of fractal parameters on shale apparent permeability evolution. These fractal parameters include nanopore radius fractal dimension, tortuosity fractal dimension and maximum nanopore radius. The analysis indicates that the variation of fractal parameters in iOM results in the change of apparent permeability of both OM and iOM, while the variation of fractal parameters in OM only leads to the change of apparent permeability of OM. The modeling results offer new insights into gas flow behaviors in shale matrix.