In this work, we develop a novel run-to-run-based model predictive controller (BR-based MPC) for a batch crystallization process with process drift and inherent variation in solubility and crystal growth rates. In order to achieve the production of crystals with desired product qualities, a conventional MPC system with nominal process model paraineters is initially applied to a batch protein crystallization process. However, the mismatch between the process model and the actual process dynamic behavior because of the process drift and variability becomes severe as batch runs are repeated. To deal with this problem of batch-to-batch variability, after each batch is over, the post-batch crystal attribute measurements, including average crystal shape and size and the number of crystals, are used to estimate off-line the drift of the process model (used in the MPC) parameters from nominal values via a multivariable optimization problem. Along with the adapted controller model parameters, the exponentially weighted moving average (EWMA) scheme is used to deal with the remaining offset in the crystal shape values and thereby to compute a set of optimal jacket temperatures. Furthermore, the crystal growth in the batch crystallization process is modeled through kinetic Monte Carlo simulations, which are then used to demonstrate the capability of the proposed R2R-based MPC scheme in suppressing the inherent variation and process drift in solubility and crystal growth rates. It is demonstrated that the production of crystals with a desired shape distribution is successfully achieved after three batch runs through the use of the proposed R2R-based MPC, while it takes 24 batch runs for the system with the EWMA-type constant supersaturation control to achieve the same objective.