There are two approaches to this question.
One of them, is separate the small scale features from the large ones, in both images, and then merge the small scales of the unsaturated version with the larse scales of the saturated one. This method is not easy to apply, so if you haven't tried this before with a single image (to perform wavelets enhancements), you should avoid this from the moment. By the way, I have used this method to combine comet exposures, using the stars as references by one hand, and the comet, by the other. So, I ended with "punctual" (point-like) stars and comet core.
The second approach is to use a propper mask. The most simple star mask (but also the less efficient) is a inverted luminance. A more refined mask may be attained with the A Trous Wavelets process. Just disable all the scales above certain scale (trial and error), for example 16px. Then, set the deringing to the maximum to all the remaining scales. This way you end with small scale features only. Some of the larger star may be appear a bit opaque. You may fix it a bit using a "large scales" image (apply to a copy of the luminance the "inverted" wavelets process, it means, disable small scales, and enable the large scales. You'll end with a blured image) as a mask, so you may increase the pixel values with HistogramTransform selectively.
Well, after you have the mask, make sure that it affects the stars, and then use PixelMath to "import" the non-saturated image. You may try to combine that version with the current one too, with a expression like this: "Saturated*0.33+NonSaturated*0.67".
Hope this helps...
