For large scale simulations you can use subdomains, e.g. layers. As far as each domain is definitely larger than onr mean free path, it should give nice results
There's no way to define features for an heterogeneous mass density for a given material, when creating an MCNP(X) input file. The only way to have something coarsely matching your request, is to define, in your material cell(s), some subcells when you can define some little variation in material bulk density, as well as in elemental composition. Of course the geometrical size for such subcells should be larger than the (e.g. neutrons) mean free path, otherwise there's no meaning. However, in case you need a fine discretization of your geometrical region, such an approach may not be done putting all data one by one, but using instead some programs able to read the different density values and automatically generate the geometrical part of your MCNP(X) input file. As an example, when estimating the flux (dose rate) map distribution for neutrons/gammas in some humans' inner organs (e.g. brain), we can start from PET/CT scans of different slices taken at different cutting levels. There are some programs able to read the raw data of PET/CT images and transform the physical information included into them, in MCNPX input file.