Most metals are crystalline materials that can undergo significant plastic (permanent) deformation when subjected to applied loading. Plastic deformation is usually accompanied by an increase in the flow stress of the material. This phenomenon is called strain hardening and is of vital importance in many engineering applications, including aerospace, automotive, and power generation industries. Developing accurate material models to predict the plastic response and hardening behavior of metals during deformation is a prerequisite to the engineering design processes, which requires a physical understanding of the underlying deformation mechanisms. In single crystals, plastic deformation of the crystal is governed by the evolution of dislocations---line defects inside the crystalline materials which marks the boundary between the slipped and unslipped regions---moving and interacting in response to the applied loading. Dislocation dynamics (DD) simulations, which track the time-space trajectories of individual dislocation lines, provide a promising tool to establish a physical link between the dislocation microstructure evolution and the strain hardening phenomenon. However, the high computational cost of DD simulations renders the accessible length and time scales to well below those which are relevant to most engineering applications. Due to this challenge, instead of directly using DD simulations for engineering applications, we have utilized DD simulations to delineate how constitutive relations of crystal plasticity (CP) can be constructed for FCC copper, based on coarse-graining of high-throughput DD simulations. This thesis consists of three main components, and we show how they fit together into a complete, physical model like three pieces of a puzzle. The first piece is a massive DD simulation database that we were able to generate thanks to recent computational advances in DD, including the subcycling time-integration algorithm and its implementation on Graphics Processing Units (GPUs). By systematically coarse-graining the database we present a strain hardening model which consists of two components: 1) a dislocation multiplication model, which accounts for slip-free multiplication, and 2) an exponential flow-rule connecting slip system shear rate to the resolved shear stress through an exponential function. These components can be thought of as the second and third puzzle pieces. By analyzing the data, it was discovered that dislocation multiplication frequently occurs on slip systems which experience zero applied shear stress (i.e., zero Schmid factor) and have a plastic strain rate of zero; we termed such multiplication slip-free multiplication and it serves as the second puzzle piece. This finding questions the assumption of the existing phenomenological expression that multiplication is proportional to the shear rate. We propose to add a correction term to the generalized Kocks-Mecking expression to account for slip-free multiplication, whose mechanistic explanation is provided. A major finding of this thesis is that DD results suggest an exponential flow-rule, in contrast to the commonly used power-law flow-rule, even in the cases where thermal fluctuations are not present. The exponential flow-rule is the third piece in the puzzle of the presented strain hardening model. We demonstrate that the observed exponential flow-rule, despite the common notion that thermal fluctuations are the responsible mechanism, can be explained by statistical properties of the dislocation links. Hence, by statistically analyzing the number density and plastic activity of links in terms of their length, we formulate a physically justified link length based flow rule which can numerically capture the exponential dependence of shear rate on shear stress. The proposed link length based flow-rule has two key components: 1) the number density of links on each slip system, which was observed to follow the sum of two exponentials distribution, and 2) an average velocity of links as a function of resolved shear stress and link length, whose fitting coefficients are independent of the loading orientation. The exponential dependence of on resolved shear stress is traced to the spatial fluctuation of the internal stress field, which can be approximated by a Laplace distribution. The proposed average velocity function incorporates the Laplace distribution in its form. This thesis shows that discrete dislocation dynamics simulations can be used to inform higher length scale models of non-phenomenological constitutive relations. The presented model captures the strain hardening as a result of slip system interactions in FCC single crystals. It works as an example for developing similar coarse-grained models based on DDD which includes additional strain hardening mechanisms such as cross-slip, or precipitate hardening. We hope that the present thesis motivates more researchers to use DDD simulations for constructing constitutive relations